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ABSTRACT  This  manual  develops  the  theory  of  the  Conditional 
Ocean  Wave  Simulation  Model  (SIMBAT).  SIMBAT  is  a  high  perform¬ 
ance  FORTRAN  F77  based  computer  model  that  uses  the  Fast  Fourier 
Transform  (FFT)  to  produce  ocean  wave  properties  for  ocean  engineering 
applications.  The  model  assumes  the  wave  properties  form  a  Gaussian 
stochastic  process.  SIMBAT  may  be  used  to  perform  a  conditional  or 
unconditional  simulation.  A  conditional  simulation  uses  a  measured  or 
existing  input  time  series  and  forces  the  simulated  wave  properties  to 
adhere  to  the  input  time  series  but  also  follow  the  laws  of  multivariate 
normal  probability.  An  unconditional  simulation  uses  a  measured  or  cre¬ 
ated  ocean  wave  spectra  to  randomly  simulate  the  ocean  wave  properties. 
Program  features  include:  creation  of  directional  ocean  wave  spectra, 
water  particle  kinematic  stretching,  conditioning  time  series  using  an  input 
time  series  less  than  or  equal  to  the  simulation  length,  use  of  Legendre 
orthogonal  polynomials  for  post-creation  of  wave  properties  by  program 
CKPOLY,  and  error  checking. 
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EXECUTIVE  SUMMARY 


This  report  provides  a  theoretical  description  of  the  Conditional  Ocean  Wave  Simulation 
Model  called  SIMBAT.  The  mathematical  foundations,  probability  theory,  and  numerical 
methods  employed  are  explained.  This  report  provides  the  user  with  the  necessary  theoretical 
information  to  utilize  SIMBAT  in  an  educated  manner.  Along  with  this  document,  a 
supplementary  document,  TN-1838  SIMBAT  User’s  Manual  (Borgman,  et  al.,  1991),  which 
provides  information  for  use  of  the  software,  is  available  from  the  Naval  Civil  Engineering 
Laboratory  (NCEL).  SIMBAT  was  developed  by  Dr.  L.  E.  Borgman  at  the  University  of 
Wyoming  for  NCEL 

The  primary  objective  of  the  SIMBAT  development  was  to  provide  an  efficient  means  to 
assist  in  analyzing  the  dynamic  motions  of  deepwater  semisubmersible  platforms  to  be  used  for 
Offshore  Tactical  Aircrew  Combat  Training  Facilities  (Shields,  et  al.,  1987).  This  development 
was  funded  by  the  Office  of  Naval  Technology  (ONT)  under  the  Navy  Exploratory  Development 
Program. 

SIMBAT  development  focused  on  the  simulation  of  water  wave  properties  to  expedite 
subsequent  calculation  of  wave  loading.  For  example,  if  the  Morison  equation  (Sarpkaya  and 
Isaacson,  1981)  is  used  to  compute  the  wave  loads  on  a  slender  member,  the  water  particle 
velocities  and  accelerations  will  be  required  as  input  for  each  time  step  and  at  each  load  point 
across  the  structure.  SIMBAT  produces  this  data  by  use  of  fast  frequency  domain  methods  and 
by  approximation  with  Legendre  orthogonal  polynomials  (Hochstrasser,  1964). 

SIMBAT  provides  the  option  to  create  either  conditionally  or  unconditionally  simulated 
water  particle  kinematics.  A  conditional  simulation  utilizes  a  measured  ocean  wave  spectrum  or 
measured  time  series  history  to  "condition"  the  simulation  to  create  associated  water  particle 
kinematics  that  adhere  to  the  properties  of  the  input  data  but  also  follow  the  laws  of  normal 
probability  theory.  This  is  particularly  useful  if  the  user  wants  to  impose  a  large  measured  wave 
profile  on  the  structure  and  create  the  associated  kinematics  for  that  profile  to  determine  the 
extreme  loads  on  the  structure.  The  advantage  of  this  option  is  that  a  large  wave  profile  is 
certain  to  occur  in  the  computer  simulation  where  ordinarily  a  much  longer  time  domain 
simulation  would  be  required  before  an  extreme  wave  profile  would  be  realizable.  Thus,  the 
design  engineer  may  utilize  SIMBAT  to  produce  the  appropriate  water  particle  kinematics  for 
extreme  wave  loading  in  a  reasonable  amount  of  computer  time. 

The  unconditional  simulation  will  produce  wave  properties  in  accordance  with  an  input 
ocean  wave  spectrum  and  that  follow  a  multivariate  normal  probability  law.  The  wave  properties 
are  "unconditionally"  simulated  randomly  by  an  input  seed  number  specified  by  the  user. 

Current  procedures  available  for  simulating  ocean  wave  properties  are  based  on  time 
domain  superposition  and  filtering  of  Gaussian  white  noise.  Substantial  savings  in  computer 
simulation  time  and  expenses  can  be  realized  if  the  wave  properties  are  simulated  in  the 
frequency  domain  as  is  done  in  SIMBAT. 

One  of  the  major  features  of  the  simulation  model  is  its  ability  to  produce  wave  properties 
over  a  large  three-dimensional  spatial  region  using  a  Legendre  polynomial  fit.  The  wave  load 
points  on  the  offshore  structure  are  calculated  at  each  time  step  in  an  expedient  manner  as 
opposed  to  the  standard  industry  method  of  computing  the  wave  properties  in  the  time  domain 
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at  each  load  point.  Orders  of  magnitude  of  computer  time  may  be  saved  as  a  result  of  this 
method.  If  SIMBAT  is  used  for  this  application,  the  design  engineer  would  merge  the  SIMBAT 
post  processor,  CKPOLY,  with  a  structural  analysis  package  to  read  the  water  particle  kinematics 
from  the  SIMBAT  output  files. 

SIMBAT  can  be  used  for  the  following  applications: 

(a)  Determination  of  ocean  wave  properties  for  computing  wave  loading  on  offshore 
structures  that  will  be  used  in  connection  with  dynamic  response  simulation 
models. 

(b)  Wave  force  studies  in  random  waves. 

(c)  Creation  of  random  directional  seas  for  design  applications. 

In  addition  to  these  three  applications,  SIMBAT  may  also  be  used  for  the  conditional 
simulation  of  waves  in  a  model  test  basin  (e.g.,  the  creation  of  a  large  wave  train  or  wave  groups 
in  a  basin  that  would  adhere  to  the  input  wave  autospectra  and  retain  the  statistical  properties  of 
the  multivariate  normal  probability  law).  This  would  allow  the  experimentalist  to  perform  tests 
in  a  shorter  period  of  time. 
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OBJECTIVE 


This  report  provides  a  theoretical  description  of  the  Conditional  Ocean  Wave  Simulation 
Computer  Model  (SIMBAT).  This  report  supplements  the  SIMBAT  User’s  Manual  (Borgman 
et  al.,  1991),  which  provides  the  instructions  for  use  of  the  SIMBAT  software. 


INTRODUCTION 

Highly  regular,  single  direction,  single  frequency  wave  trains  are  quite  rare  in  the  ocean. 
The  predominant  wave  systems  present  are  irregular,  with  at  least  some  degree  of  short 
crestedness  (Figure  1).  The  wave  conditions  in  an  intense  storm  are  particularly  chaotic  and 
turbulent,  yet  these  are  the  conditions  most  significant  to  the  engineer  designing  a  structure  for 
emplacement  in  the  ocean. 

Several  options  are  available  to  the  designer.  Nonlinear,  unidirectional  waves  of 
permanent  form,  such  as  those  produced  by  Stokes’  theory  or  various  numerical  solutions,  can 
be  used  for  each  large  "bump"  of  sea  surface.  This  produces  a  conservative  or  maximal  estimate 
of  the  water  particle  velocities  and  accelerations,  because  each  component  of  the  wave  form 
reinforces  the  others  to  push  and  pull  together.  Alternatively,  a  linear  superposition  of  many 
independent  components,  each  with  its  own  direction,  phase,  and  amplitude,  can  be  used  to 
obtain  the  wave  kinematics. 

Both  of  these  options  are  imperfect  models  of  the  real  ocean.  The  fault  of  the  nonlinear, 
unidirectional  theory  is  that  it  imposes  an  artificial  regularity  on  its  estimates  of  the  wave 
kinematics.  The  theory  based  on  linear  superposition  allows  the  irregularity,  but  in  turn  produces 
a  linear  solution  to  wave  conditions  that  are  known  to  have  skewed  nonlinear  characteristics  of 
sharper  crests  and  flatter  troughs.  The  design  engineer  must  judge  which  approach  has  the  least 
imperfection  for  the  application  intended.  There  are  situations  where  each  is  the  most 
appropriate  compromise.  A  discussion  of  these  options  and  a  comparison  with  field  data  are 
given  by  Forristall,  Ward,  Borgman,  and  Cardone  (1978). 

More  to  the  point,  the  theory  used  for  the  wave  kinematics  must  be  evaluated  in  tandem 
with  the  wave  force  formula  and  the  values  of  the  force  coefficients  selected.  An  overestimation 
within  the  wave  kinematics  algorithm  can  be  somewhat  compensated  for  with  lower  force 
coefficients,  and  vice  versa. 

The  Morison  equation  is  generally  accepted  as  a  satisfactory  approximation  for  wave 
force,  although  other  alternatives  and  variations  have  been  investigated  from  time  to  time.  One 
common  modification  is  to  vary  the  drag  coefficient  with  time  according  to  some  locally 
evaluated  Reynolds  number. 

The  particular  combination  of  wave  theory  and  force  formula  selected  for  use  really  needs 
to  be  calibrated  against  measured  wave  data  under  real  oceanic  conditions.  Most  oil  companies 
have  gone  through  this  exercise  at  least  once  for  the  schemes  they  have  found  useful  in  design. 
The  larger  companies  proceed  with  such  evaluations  almost  continually  to  match  new  oceanic 
conditions  as  they  arise  in  the  company’s  operations.  There  is,  of  course,  some  diversity  in 
computational  procedures  from  company  to  company,  because  each  may  choose  to  adjust  the 
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overall  kinematics  and  force  calculations  in  different  places  to  calibrate  the  produced  estimates 
to  agree  with  measured  data. 

Various  reviews  of  nonlinear  wave  theory  (Dean,  1970;  Sarpkaya  and  Isaacson,  1981), 
random  wave  theory  (Borg  man,  1972;  Ochi,  1982)  and  force  formulas  (Dean  and  Borg  man, 
1986;  Sarpkaya  and  Isaacson,  1981)  have  been  published.  Dean  and  Borgman  (1986,  p.  341) 
also  give  a  good  summary  of  field  measurement  programs.  It  does  not  appear  useful  to  provide 
yet  another  summary  of  these  topics.  Instead,  various  new  results  largely  related  to  the 
estimation  of  kinematics  in  irregular  wave  fields  as  input  to  force  computations  will  be  presented. 
Computer  simulation  of  wave  kinematics  and  the  special  topic  of  conditional  simulation  will  be 
particularly  emphasized.  The  inclusion  of  the  word  "forces"  in  the  chapter  title  is  intended  to 
indicate  that  the  formulas  for  kinematics  are  directed  toward  providing  a  basis  for  wave  force 
computations.  New  public  domain  computer  software  is  just  now  becoming  available  that 
enhances  the  usefulness  of  these  techniques  in  engineering  design.  Sources  for  such  software  will 
be  briefly  reviewed. 


Figure  1 

Aerial  photograph  of  an  irregular  sea  surface  (ship  in  photo 
shows  scale)  (from  Pierson,  Neumann,  and  James,  1990). 
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BACKGROUND 


The  Naval  Civil  Engineering  Laboratory  (NCEL),  under  the  Offshore  Aircrew  Combat 
Training  System  task  of  the  Navy  Exploratory  Development  Program,  was  tasked  to  develop 
technology  necessary  to  design  and  construct  reliable  and  cost  effective  unmanned  ocean 
platforms  for  deepwater  (600  to  10,000  feet)  Tactia!  Aircrew  Combat  Training  System  (TACTS) 
applications  (Shields,  et  al.,  1987).  One  of  the  objectives  of  this  development  program  was  to 
develop,  modify,  integrate,  and  validate  computer  models  for  simulation  of  deepwater  moored 
platform  motions.  Investigations  of  various  candidate  platforms  resulted  in  the  semisubmersible 
concept  being  selected  as  the  most  suitable  platform  type  (Shields,  et  al.,  1983). 

The  SEASTAR  (PMB  Engineering,  1990)  nonlinear  structural  finite  element  analysis 
program  was  selected  to  model  the  hydrodynamic  loads  and  dynamic  response  of  the 
semisubmersible/ mooring  system.  SIMBAT  was  developed  to  meet  the  following  requirements: 

1 .  Wave  force  studies  in  random  seas. 

2.  Wave  property  input  in  dynamic  response  simulation  model,  SEASTAR. 

3.  Random  directional  seas  representation  for  design  applications. 

In  order  to  gain  an  insight  into  the  dynamic  response  of  deepwater  TACTS 
semisubmersible  platforms,  the  dynamic  response  simulation  model  SEASTAR  was  used.  The 
preferred  environmental  data  base  for  SEASTAR,  when  used  in  design  analyses,  consists  of  a 
number  of  actual  ocean  wave  time  series  corresponding  to  major  storms  (i.e.,  50,  75,  100  year 
storm)  or  other  environmental  scenarios  of  interest  (i.e.,  those  for  fatigue  analyses  or  operational 
sea  states)  at  the  sites  of  interest.  Since  these  data  are  unavailable  in  most  design  projects, 
simulated  ocean  wave  property  data  should  be  used. 

Current  procedures  available  for  simulating  ocean  wave  properties  are  based  on  time 
domain  superposition  and  filtering  of  Gaussian  white  noise.  Substantial  savings  in  simulation 
time  and  expenses  can  be  realized  if  the  discrete  Fourier  transform  of  the  desired  time  histories 
of  the  wave  properties  are  simulated  directly  in  the  frequency  domain.  The  frequency  domain 
simulation  methods  are  ten  to  a  hundred  times  faster  than  the  time  domain  simulation  techniques. 
In  addition,  the  computationally  very  rapid  Fast  Fourier  Transform  (FFT)  can  be  used  to  revert 
the  simulation  to  the  time  domain. 

All  simulation  schemes  are  based  on  the  introduction  '■‘f  random  numbers.  The 
engineer/oceanographer/mathematician  builds  in  the  wave  properties  (design  wave/wave 
grouping)  to  be  preserved.  If  one  has  an  actual  sequence  of  ocean  wave  data  and  wishes  to  study 
the  motion  response  associated  with  that  particular  sequence,  currently  available  simulation 
procedures  will  not  accommodate  these  data.  Conditional  simulations  are  necessary,  with  the 
randomness  being  restricted  so  as  to  produce  ocean  wave  simulation  data,  conditional  on  the 
wave  having  the  required  values.  These  required  values  may,  for  example,  be  a  series  of  wave 
heights  expected  in  a  typical  storm  or  other  environmental  scenarios  at  the  site. 

Conditional  simulation  of  ocean  wave  properties  is  a  very  recent  research  topic  with  wide 
application  to  Navy  projects.  Techniques  of  conditional  simulation  have  been  used  in  geological 
problems  and  are  now  applied  to  ocean  wave  applications  as  used  in  SIMBAT. 

During  the  development  of  SIMBAT,  the  use  of  Legendre  orthogonal  polynomials 
(Hochstrasser,  1964)  was  added  to  SIMBAT  to  create  three-dimensional  wave  properties  over 
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a  large  spatial  region.  The  use  of  these  polynomials  along  with  the  frequency  domain  method 
of  creating  the  wave  properties  was  determined  to  provide  a  computationally  much  faster 
algorithm  than  computing  the  wave  properties  at  all  grid  points  over  all  time.  This  application 
utilizing  Legendre  polynomials  is  now  considered  one  of  the  major  beneficial  aspects  of 
SIMBAT,  especially  for  compliant  systems. 


OVERVIEW  OF  SIMBAT  MODEL  AND  SOFTWARE 

SIMBAT  is  a  computer  model  that  calculates  ocean  wave  properties,  such  as  water 
particle  kinematics  for  use  in  various  ocean  engineering  applications.  The  model  assumes  the 
wave  properties  form  a  Gaussian  stochastic  process.  This  means  that  all  wave  properties  at  any 
selected  set  of  times  and  locations  follow  a  multivariate  normal  probability  law. 

The  SIMBAT  software  can  create  the  Ochi-Hubble,  Pierson-Moskowitz,  Bretschneider, 
JONSWAP,  and  Wallops  Spectral  Models  or  read  in  the  user’s  own  spectra,  if  desired.  The 
waves  may  be  directionally  spread  by  either  wrapped  normal,  cosine-squared,  or  von  Mises 
methods.  There  are  options  to  stretch  wave  properties  above  mean  water  level  using 
Rodenbusch-Forristall  delta  stretch,  Reid- Wheeler  stretch,  truncation  extrapolation,  functional 
extrapolation,  linear  extrapolation,  and  gamma  extrapolation. 

SIMBAT  has  the  option  to  perform  a  single  si  mu'  ition  of  wave  properties  at  one  or  more 
location:  in  space  using  frequency  domain  methods  directly  or  using  the  application  of  Legendre 
orthogonal  polynomials.  If  the  polynomials  are  used,  a  separate  program,  CKPOLY,  is  used  to 
create  the  wave  properties  from  SIMBAT  output.  CKPOLY  is  provided  with  SIMBAT  but  is 
separate  to  allow  the  user  to  embed  the  program  in  their  structural  analysis  software. 

For  conditional  simulation,  the  user  can  provide  SIMBAT  with  a  measured  large  wave 
profile,  wave  groups,  or  any  other  wave  properties  that  are  presentative  of  the  environmental 
area  of  interest.  SIMBAT  will  constrain  the  simulated  wave  properties  to  follow  the  measured 
data  while  also  maintaining  the  statistical  laws  appropriate  for  ocean  waves. 

The  latest  version  of  SIMBAT,  Release  3.0,  now  contains  accuracy  verification 
algorithms  to  compare  kinematics  from  the  different  simulation  methods  in  SIMBAT. 

SIMBAT  has  the  option  to  write  ASCII  output  data  files  as  it  proceeds  through  the 
simulation  so  that  the  user  may  verify  the  integrity  of  the  data.  Therefore,  if  SIMBAT  is 
executed  in  a  multitasking,  multiwindow  environment,  the  user  can  confirm  the  simulation  is 
working  properly  on  his/her  computer  throughout  each  level  of  the  simulation. 

This  study  was  motivated  by  the  need  for  computer-efficient  calculation  procedures  and 
computer  codes  for  the  preparation  of  accurate  ocean  wave  kinematics  and  properties  for 
multiple-frequency,  multiple-direction  seas.  The  immediate  application  was  the  design  of 
moored,  floating  instrument  towers  for  Navy  gunnery  ranges.  However,  the  extensive  theory 
and  efficient  computational  algorithms  should  prove  useful  in  many  other  ocean  engineering 
applications. 
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THEORETICAL  FORMULATIONS  OF  THE  SBMBAT  MODEL 


Coordinate  System  Specifications 

The  ocean  wave  kinematics  will  be  referenced  to  a  general  horizontal  coordinate  system. 
All  horizontal  coordinate  axes  are  established  within  navigation  headings  measured  clockwise 
from  true  north: 

0X  =»  direction  of  positive  x  axis 

8  =  direction  of  positive  y  axis  (1) 


0-8  **  90° 

x  y 


Let  the  vertical  axis  z  be  zero  at  mean  water  level  and  positive  downward. 

The  direction  of  travel  of  a  wave  is  8  in  navigation  heading.  The  wave  is  traveling 
toward  direction  8  if  p  =  1  and  is  coming  from  direction  8  if  0  —  -1. 

Basic  Wave  Properties 

Eight  wave  oroperties  a^  of  interest.  In  terms  of  «eal  functions,  these  are: 

1.  The  water  level  elevation, 

ri(x,y,t)  =  a  cos  {  P0k[x  cos  (8  -0X)  +y  cos(8  -6y)]  -  2  it  ft  -  4>)  (2) 

2.  The  components  of  water  particle  velocity, 


Vx(x,y,z,t) 

-  cosh[k(d-z)] 

p0  cos(0  -  8X) 

Vy(x,y,z,t) 

sinh  (k  d) 

p0  cos(0  —  8y) 

Vz(x,y,z,t)  =  a(2af) 


*  cos[p0k{x  cos  (8  -0X)  *y  cos(8  -0,,)}  -  2  a  ft  -  <j>] 
sinh[k(d-z)] 


sinh(kd) 

*  sin  [  P0 k  { x  cos  (8  -  0X)  +  y  cos  (0  -  0y)}  -  2  a  f  t  -  4>]  ^4^ 
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3. 


The  components  of  water  particle  acceleration, 


a,(x,y,z,t) 

[ay(x,y,z,t) 


a(2*f)!  cogh[k(d-z)] 
sinh(kd) 


p0  cos  (0  -  0,) 
P0  co6  (0  -  0y ) 


*  sin[P0k|x  cos(0  -0x)  +  y  cos(0  -0y))  -  2  it  ft  -  $1 


(5) 


ax(x,y,z,t) 


-  a  (2  n  f )J  sMi1I[(J-2» 
sinh(kd) 


*  cos  [  P0  k  { x  cos  (0  -  0x)  +  y  cos  (0  -  0y)  1  -  2  n  f  t  -  <J>]  (6) 


4.  The  water  pressure  anomaly  (plus  and  minus  about  hydrostatic  pressure)  divided  by  pg, 

pfou.M)  =  a  cosh[k(d  -  z)] 
pg  cosh(kd) 

*  cos  [  P0  k  { x  cos  (0  -  0x)  +  y  cos  (0  -  0y) }  -  2  n  f  t  -  4>]  (7) 


In  these  formulas: 

a  =  wave  amplitude 
f  =  wave  frequency 
d  =  water  depth 

k  =  wave  number  =  27r/wave  length 
0  =  wave  phase 
p  =  water  density 
g  =  acceleration  due  to  gravity 

and  it  is  assumed  that  wave  number  and  wave  frequency  are  related  by  the  formula: 

(2  it  f)2  =  gktanh(kd)  (8) 


Wave  Properties  in  Complex  Form 

Through  the  use  of  the  complex  form  and  cos  a  and  sin  a,  where 


cost  =  e*p(i«>  -  exp  (-  i  a) 

2 


(9) 
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sina  _  exp-flq)  -  exp(-i«) 

2i 

All  of  the  wave  properties  listed  above  can  be  expressed  in  the  form: 


b(x,y,z,t) 


G(z) T(f)  H(0)  exp[-  i  p0 k  { x  coe(6  -  0,) 


y  cos(0  -  0y)J  ]  exp(i  2  n  f  t) 


for  positive  f.  The  original  real-valued  wave  time  property  equals  b(f)  ■+•  b(f)*,  where 

b*(f)  =  complex  conjugate  of  b(f) 

The  functions  G,  T,  and  H  for  each  wave  property  are: 

1.  Water  level  elevation, 

G(z)  =  T(f)  =  H(0)  =  1.0 


2.  Velocity, 


cosh[k(d-z)] 
sinh  (kd) 

sinh  [k(d  -  z)] 
sinh(kd) 


{e~kl  +e*k(2d'2)} 


{1  -  e'Zkd} 

{e  kr  -  e  k(2d  l)} 
{1  -  e~2kd} 


- -  ,  for  V  and  V 

-2kd>  *  *  y 


,  for  V, 


2  n  f  , 

for  V, 

and  Vv 

A 

y 

(15) 

2  n  f  i  , 

for  V, 

Pocos(0 

forV, 

Pocos(0 

-ey). 

for  Vy 

(16) 

1.0, 

for  V, 
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3.  Acceleration, 


G(z) 


cosh  [k(d  -  z)] 
sinh(kd) 


< 


{e'kl  -t-e  tod  ,)} 
{1  -  e'2kd) 


sinh[k(d  -z)]  _  {e~k*  - 

sinh(kd)  |l-e'2kd) 


for  az  and  ay 
for  a. 


(17) 


4. 


Pressure  anomaly, 

G(z) 


T(f)  = 

(2  n  f  )2  i  , 

~  (2  it  f  J2  , 

for  az  and  ay 

,  for  az 

08) 

P0  cos (0 

-  0Z)  ,  for  az 

H(0)  = 

P0  cos  (0 

-  0y)  ,  for  ay 

(19) 

1.0  , 

for  az 

cosh  [k(d  -  z)] 

(e-k*  +e-k(2d-D| 

(20) 

cosh(kd) 

(1  +  e~2kd} 

II 

© 

H 

1.0 

(21) 

H(0)  - 

1.0 

(22) 

In  the  above  definitions,  the  hyperbolic  function  terms  involving  the  vertical  coordinate 
z  and  water  depth  d  have  been  expressed  in  algebraically  equivalent  forms,  which  makes  it  easy 
to  extend  the  formulas  to  very  large  water  depths.  These  forms  also  have  advantages  later  in  the 
sections  on  the  treatment  of  compliant  structures. 

Irregular  Waves  by  Superposition 

An  irregular  wave  train  can  be  obtained  by  the  summation  of  many  different  simple  wave 
forms.  For  example,  the  general  formulas  in  Equation  11  can  be  summed  to  obtain,  with  some 
obvious  extensions  of  notation, 
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(23) 


b(t*,y,z)  -EE  «.,e‘*-'G(f„,z)T(f„)H(8) 

m-  M  j-1 

*  expI-iPok^lx  cos (0J  - 6,)  +  y  cos^  -0y)}]exp(i2nfmt) 


The  summation  extends  over  an  arbitrarily  selected  list  of  frequencies  {fm;  m  =  1,2,3,...,M}, 
and  a  corresponding  list  of  directions  (0j;  j  =  1,2,3,...,J}.  By  convention, 

f0  -  0.  f-  -  -fm  <24) 


and  the  terms  in  the  summation  for  negative  frequencies  are  taken  to  be  the  complex  conjugates 
of  those  for  the  corresponding  positive  frequencies.  This  forces  the  wave  property  time  series 
to  be  real-valued. 

Irregular  Waves  in  Frequency  Domain 
Let 


B<f„;w)  =  E  ^,e,*"G(f„,z)T(f„)H(8|) 

j-i  (25) 

*  expI-iPokJx  cos^-OJ  +  y  cos^-Oy)}] 

Equation  23  may  be  expressed  as  a  discrete  Fourier  transform, 

M 

b(t;x,y,z)  =  £  B(fm;x,y,z)exp(i2  n  fmt)  (26) 

m«-M 


This  form  is  particularly  convenient  since  it  permits,  with  some  slight  modifications,  the  high¬ 
speed  computation  of  wave  property  time  series  with  the  Fast  Fourier  Transform  (FFT) 
algorithm. 

Fast  Fourier  Transform  (FFT)  Algorithm 

The  FFT  algorithm  provides  a  very  efficient  and  rapid  procedure  for  calculating  the  two 
formulas  (Blahut,  1985): 
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w 


(27) 


N-l 


“  E  w-exp(-i2nmn/N) 

N  q.o 


N-l 


w„  -  E  Waexp(i2nmn/N) 

m«0 


(28) 


for  m  =  0,1,2,. ..,N-1  and  n  =  0,1,2,. These  are  exact  discrete  transformations  for 
computing  one  of  the  two  sequences,  {Wm;  m  =  0,1,2,...,N-1}  and  {wn;  n  =  0,1,2,...,N-1}, 
from  the  other. 

Some  properties  of  the  FFT  formulas  deserve  mention.  Because 

exp(±  i2  «  mn/N)  —  cos (2  it  mn/N)  ±  isin(2  n  mn/N)  (29) 


both  summations  are  formally  periodic,  with  period  N.  This  means  that  both  sequences  {Wm} 
and  {Wn}  are  forced  by  their  mathematics  to  have  this  periodicity.  Also,  both  sequences  are 
complex-valued,  in  general.  However,  in  the  applications  here,  the  wn  are  real-valued  and  the 
Wm  are  complex-valued.  A  necessary  and  sufficient  condition  for  the  wn  to  be  real-valued  is 
that: 


W 


m 


Wn.* 


(30) 


for  0  <  m  <  N/2. 

Wave  Properties  in  the  FFT  Format 

The  conversion  c  ”  Equation  23  to  the  FFT  format  is  achieved  by  discretizing  time  and 
frequency, 


t  =  nAt,  n  =  0,1, 2,. ...N-l 

fm  =  mAf,  m  =  0,1,2,..., N/2 


(31) 


where  the  time  and  frequency  increments  are  taken  to  satisfy 

(At)(Af)  -  1/N  (32> 

With  this  constraint,  the  argument  in  the  complex  exponential  in  Equation  25  simplifies  to: 

2  x  fm  t  -  2  x  m n/N  (33) 

In  addition,  the  periodicity  implicit  in  the  FFT  formulas  leads  naturally  to  the  translation  of  the 
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negative  frequency  values  in  Equation  26  to  the  positive  frequencies: 

B({N-m}  Af;x,y,z)  -  B(-m Af;x,y,z)  C*4) 

The  complex-valued,  skew  conjugate  symmetry  imposed  on  Equation  23  then  becomes  the  FFT 
symmetry  in  Equation  30  that  forces  the  transform  to  be  real-valued. 

With  these  modifications,  the  general  wave  property  formulas  in  Equation  26  become: 

N-l 

b(nAt;x,y,z)  =  £  B(m  Af;x,y,z)exp(i2  n  mn/N)  (35) 

m»0 

The  values  of  N  and  Af  should  be  selected  so  that  all  the  energetic  frequencies  of  importance  are 
enclosed  within  the  bounds: 

0  <  mBAf  <  f  <  mLAf  <<  1/(2  At)  (36) 


The  interval  (m^m^  ordinarily  extends  over,  perhaps,  150  Af  increments  out  of  N  =  2048 
terms  in  the  sequence.  The  other  FFT  coefficients  can  be  set  to  zero,  except  for  those 
determined  by  complex  conjugation  in  symmetry  about  m  =  N/2. 

In  summary,  a  frequency  domain  computation  of  the  wave  property  time  series  consists 
of  the  following  steps: 

1.  B(mAf;  x,  y,  z)  is  computed  from  Equation  25  for  mB  ^  m  £  mL. 

2.  By  conjugate  symmetry,  for  mB  £  m  s  mL, 

B({N  -m)  Af;x,y,z)  *  B*(m  Af;x,y,z)  <37) 

3.  The  rest  of  the  coefficients  are  set  to  zero. 

4.  The  sequence  of  B(mAf;  x,  y,  z)  values  are  transformed  with  the  Fast  Fourier  Transform 

algorithm  by: 


N-l 

b(nAt;x,y,z)  =  £  B(m  Af;x,y,z)exp(i2  n  mn/N)  (38) 

m-0 

This  total  process  yields  the  full  time  series  for  n  =  0,1, 2,..., N-l. 
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THE  COMPLEX- VALUED  AMPLITUDE  MATRIX 

The  irregular  train  and  all  its  associated  linear  properties,  such  as  its  kinematics,  are 
specified  completely,  for  the  FFT  approach,  by: 

A.,  -  •.,«“«  <39> 

for  mB  <  m  <  mL  and  j  =  1,2,...,J.  The  amplitudes  Amj  can  be  arranged  into  a  complex  matrix 
with  rows  specifying  direction  and  columns  designating  frequencies.  The  other  functions  in 
Equation  23  are  transfer  functions  that  convert  the  amplitudes  into  the  various  wave  properties. 

Any  procedure,  either  deterministic  or  random,  that  generates  the  complex  amplitude 
matrix  fully  determines  all  the  wave  property  time  series.  Much  of  the  following  discussion  will 
be  concerned  with  various  choices  and  their  consequences  related  to  the  probabilistic  selection 
of  values  for  the  complex  amplitudes. 

Characterization  of  Irregular  Waves 

The  directional  spectral  density  S(f,6)  is  the  function  most  commonly  used  to  characterize 
the  mean-square  oscillation  of  the  wave  components  at  different  frequencies  and  directions. 
More  specifically,  define  the  two-sided  spectrum  as: 

-  ai,/2 

=  mean-square  oscillation  of  a  cosine 

wave  with  amplitude  amj  (40) 


The  mean-square  oscillation  of  the  sea  surface  in  the  irregular  wave  train,  if  all  components 
behave  independently,  becomes: 

E  E  -2/  /  stf.ejdedf  <«) 

n»-l  J»1  2  0  0 

In  the  case  where  the  sea  surface  is  taken  as  a  random  process,  this  becomes: 

"  2* 

var(q)  =  2  f  f  S(f,0)d0df  (42) 

o  o 

A  function  of  angle  0  at  each  frequency,  called  the  spreading  function,  is  defined  as: 
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D(0,f)  -  S(f,6)/S(f) 


(43) 


where 


S(f)  =  f  S(f,0)de 


The  relation  in  Equation  43  can  also  be  written  as: 

S(f,0)  =  S(f)D(0;f) 


Many  formulas  have  been  proposed  and  studied  as  models  for  S(f)  and  the  spreading 
function  (Borgman,  1979;  Ochi,  1982).  The  most  effort  has  been  spent  seeking  the  best 
characterization  for  fully  arisen  seas,  in  order  to  improve  wave  forecasts  and  hindcasts.  The 
emphasis  for  structural  design  is  somewhat  different,  where  it  is  desirable  to  be  able  to  model 
a  wide  variety  of  sea  conditions,  from  very  narrow-banded  to  very  broad-banded,  and  from 
unidirectional  conditions  out  to  broadly  spread  seas.  Of  particular  importance  are  seas  that  are 
multimodal  in  direction  and/or  frequency. 

All  of  the  spreading  functions  commonly  used  (the  cosine-squared,  the  von  Mises,  and 
the  wrapped  normal  formulas)  have  almost  exactly  the  same  shape.  That  is,  parameters  can  be 
selected  that  make  the  formulas  very  nearly  lie  on  top  of  each  other.  All  of  these  are  unimodal 
and  symmetric  about  the  principal  direction. 

There  are  several  very  good  discussions  of  these  various  choices  of  formulas  (Sarpkaya 
and  Isaacson,  1981,  pp.  504-520;  Ochi,  1982,  pp.  308-346;  and  Muga,  1984,  pp.  163-175). 
Rather  than  repeating  still  another  summary  here,  an  often  overlooked  specific  choice  that  has 
many  advantages  in  engineering  studies  will  be  introduced.  This  choice  consists  of  the  Ochi- 
Hubble  spectrum,  as  combined  with  the  wrapped  normal  spreading  function.  Superpositions  of 
such  spectral  models  can  be  used  to  obtain  multimodal  conditions. 

Ocean  Spectral  Models 

Ochi-Hubble  Spectra  Formula.  The  formulation  developed  by  Ochi  and  Hubble  (1976) 
was  expressed  as  a  one-sided  spectrum  with  radian  frequency.  For  consistency  with  the  present 
treatment,  this  is  changed  to  a  two-sided,  cycles-per-second  formula, 

2[(4 A  ♦  Dfp/4]1  oIexp[-(4A.  ^  l)(f0/f)4/4]  (46) 

r(X)  |f|4A+1 

where  r(X)  is  the  complete  gamma  function.  It  is  straightforward  to  show  by  calculus  that  S(f) 
integrates  over  (-<»,  »)  to  obtain  a2.  The  change  of  variable, 
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y  -  (4X+l)(f0/f)4/4 


(47) 


is  made  within  the  integral  to  reduce  the  expression  to: 


variance  =  2  J  S(f)df  =  o^TfX)]'1  J  e'yyA",dy  (48) 

o  o 

If  S(f)  is  differentiated  with  respect  to  f  and  set  equal  to  zero,  to  determine  the  value  at  which 
S(f)  is  a  maximum,  the  result  is  f  =  f0.  The  maximum  value  for  the  spectrum,  S0,  is  then: 

$0  -  2.L(.t^.i 1 )/ 4.)*.  q2  exp[  - (4  X  + 1)/4]  (49) 

0  r(X)f0 


From  this, 


p  =  Vo  7[(4  X  +  l)/4]*  exp[-(4  A,  +  1 )/ 4]  = 
o2  T(X) 


(50) 


where  G(X)  is  the  function  of  X  and  P  the  parameter  so  defined. 

The  effective  width  6  of  a  spectral  density  will  be  defined  as  the  width  of  a  square  pulse 
with  the  same  height,  S0,  and  area,  a2 12,  as  the  spectral  density  from  (0,  «>).  That  is, 


6S0  =  — 
0  2 


(51) 


6  = 


2Sn 


(52) 


Another  expression  for  the  parameter  P  is: 


P  =  _L  =  G(X) 
25 


(53) 


This  conveniently  relates  X  to  the  effective  width  and  peak  frequency  of  the  spectrum  formulas. 

If  X  =  1,  the  Ochi-Hubble  spectrum  reduces  to  the  usual  Pierson-Moskowitz- 
Bretschneider  formulas.  If  X  >  1,  say  100  or  so,  the  spectrum  is  very  narrow-banded.  If  X  < 
1,  the  spectrum  is  broader  than  the  usual  spectral  models.  Table  1  shows  P  versus  X.  It  is  easy 
to  generate  more  elaborate  tables  from  Equation  50  using  a  computer,  or  to  solve  by  Newton- 
Raphson  iteration  (Press,  et  al.,  1986,  pp.  254-259)  for  the  value  of  X  required  to  achieve  a  give 
P  value. 
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Table  1 

Lambda  Versus  Spectral  Parameters  in  the  Ochi-Hubble  Formula 


Lambda/ 

Multiplier 

Multiplier  | 

0.01 

0.1 

1.0 

10.0 

100.0 

1.0 

1 

/fea 

0.716 

2.49 

7.97 

1.1 

JH 

0.758 

2.62 

8.36 

1.2 

0.0183 

0.156 

0.798 

2.74 

8.73 

1.3 

0.0198 

0.167 

0.836 

2.85 

9.09 

1.4 

0.0213 

0.177 

0.873 

2.96 

9.43 

1.5 

0.0228 

0.188 

0.908 

3.07 

9.77  f 

1.6 

0.0242 

0.198 

0.942 

3.17 

10.09  ) 

1.7 

0.0257 

'  - | 

0.975 

3.27 

10.40 

1.8 

0.0272 

0.006 

3.336 

10.70 

1.9 

0.0286 

0.227 

1.037 

3.46 

10.99 

2.0 

0.0301 

0.237 

1.067 

3.55 

11.28 

2.1 

0.0315 

0.246 

1.096 

3.64 

11.56 

2.2 

0.0330 

0.255 

1.125 

3.72 

11.83 

2.3 

0.0344 

0.264 

1.153 

3.81 

12.09 

2.4 

0.0359 

0.273 

1.180 

3.89 

12.35 

2.5 

0.0373 

0.281 

1.206 

3.97 

12.61 

2.6 

0.0387 

0.290 

1.232 

4.05 

12.86 

2.7 

0.0401 

0.298 

1.258 

4.13 

13.10 

2.8 

0.0416 

0.306 

1.283 

4.20 

13.35 

2.9 

0.0430 

0.314 

1.307 

4.28 

13.58 

3.0 

0.0444 

0.322 

1.331 

4.35 

13.81 

3.2 

0.0472 

0.338 

1.378 

4.50 

14.27  | 

3.4 

0.0500 

0.353 

1.423 

4.64 

14.71 

3.6 

0.0527 

0.368 

1.467 

4.77 

15.13 

3.8 

0.0555 

0.382 

1.510 

4.90 

15.55 

4.0 

0.0582 

0.396 

1.551 

5.03 

15.% 

4.2 

0.0609 

0.410 

1.592 

5.16 

16.35 

4.4 

0.0636 

0.423 

1.631 

5.28 

16.73 

4.6 

0.0663 

0.436 

1.670 

5.40 

17.10 

4.8 

0.0690 

0.449 

1.707 

5.51 

17.48 

5.0 

0.0717 

0.462 

1.744 

5.63 

17.83  j 

5.5 

0.0782 

0.492 

1.833 

5.91 

18.70 

6.0 

0.0847 

0.521 

1.918 

6.17 

19.54 

6.5 

0.0911 

0.548 

1.999 

6.42 

20.34 

7.0 

0.0974 

0.575 

2.077 

6.66 

21.11 

7.5 

0.1036 

0.600 

2.152 

6.90 

21.86 

8.0 

0.1097 

0.625 

2.225 

7.13 

22.57  ! 

8.5 

0.1157 

0.649 

2.295 

7.35 

23.27  ! 

9.0 

0.1217 

0.672 

2.364 

7.56 

23.95 

9.5 

0.1276 

0.694 

2.430 

7.77 

24.60 

9.9 

0.1322 

0.712 

2.482 

7.93 

25.09 

"Tabled  value  =  peak  frequency*spectral  peak/variance  =  peak  frequency/(2*effective  width). 
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Figure  2 

Decomposition  of  a  bimodal  spectral  density. 


As  suggested  by  Ochi  and  Hubble  (1976)  the  formula  is  particularly  useful  for  the 
approximation  of  multimodal  spectra.  The  formula  also  appears  to  have  a  much  wider  range  of 
applicability  than  was  investigated  by  the  authors.  Consider  an  example,  shown  in  Figure  2, 
where  the  variance  of  the  total  spectrum  is  taken  to  be  a2  =  5m2.  Suppose  2/3  of  the  variance 
is  ascribed  to  the  main  mode.  Then, 


f01  =0.1  Hz,  S01  =  33.0,  a]  -  ^ 

*02  =  0.2  Hz,  $02  =  8.0,  o\  - 


(54) 


and 

P,  =  0.990,  P2  =  0.960 


(55) 


For  these  values,  the  k  parameters  are: 

A.,  =  1.748,  X2  =  1.655 


(56) 


The  total  spectrum  is  then  modeled  with  a  combination  of  these  two  modes. 

In  practice,  modes  are  added  to  the  model,  the  residuals  examined,  and  further  modes 
introduced  until  an  adequate  approximation  is  obtained.  If  an  f5  spectral  slope  is  desired  for 
high  frequencies,  a  mode  with  k  =  1  can  be  added  to  the  right-hand  tail  of  the  spectrum,  while 
all  the  other  modes  at  lower  frequencies  are  kept  much  narrower. 
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Pierson-Moskowitz  Spectrum  (Pierson  and  Moskowitz,  1964).  The  Pierson-Moskowitz 
spectrum  is  used  extensively  by  ocean  engineers  as  one  of  the  most  representative  for  waters  all 
over  the  world  (Sarpkaya  and  Isaacson,  1981;  Chakrabarti,  1987): 


S(f) 


g  g2  t-i-zstr/y4) 

(2  it)4f5 


where:  a  =  8.1  x  10'3  =  Phillips’  constant 

f0  =  peak  frequency  (Hz) 
g  =  gravity  (ft/s2) 

Bretschneider  Spectrum  (Bretschneider,  1959).  The  Bretscheider  spectrum  is  designed 
to  ensure  the  area  mQ  under  the  spectrum  corresponds  to  H/16,  which  assumes  a  Rayleigh 
distribution  of  wave  heights  (Sarpkaya  and  Isaacson,  1981).  The  Bretscheider  spectrum  is  based 
on  the  assumption  that  the  spectrum  is  narrow-banded  and  the  individual  wave  height  and  wave 
period  follow  the  Rayleigh  distribution  (Chakrabarti,  1987): 


S(f) 


5H? 

16  f. 


'_1_J 

,<f/fo>5, 


eI(-5/4)(f/f0)-4J 


(58) 


where:  fQ  =  peak  frequency 

Hs  =  significant  height 

The  relationship,  Ts  =  0.946  T0,  where  Ts  =  significant  wave  period,  makes  the 
Bretschneider  and  Pierson-Moskowitz  (P-M)  models  equivalent  (Chakrabarti,  1987). 

JONSWAP  Spectrum  (Hasselman,  et  al.,  1973).  JONSWAP  is  a  modification  of  the 
P-M  spectrum  to  account  for  the  effect  of  fetch  restrictions  and  to  provide  for  a  much  more 
sharply  peaked  spectrum  (Sarpkaya  and  Isaacson,  1981;  Barltrop  and  Adams,  1991): 


SCO 


g  g2  e«-3/4)(f/y 4j i* 

(2  71  )4  f5 


(59) 
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where: 


i-v-tfnft1'] 


a  =  e 


o§  =  0.07  for  f  s  fo 


ob  =  0.09  for  f  >  f0 


H.  ft 

603.9 

A  8 


,2  \2  036 


(1.0  -  0.298  In  Y) 


Wallops  Spectrum  (Huang,  et  al.,  1981). 


where: 


S(f)  =  P82  eK-®/W04l 


log(y/2  n  e)2 
log  2 


- 

(2  n  e)2m(i/4)(tn~l) 

4(l/4)(m-5) 


4)(m-l)  r  1  | 

*  |r[(l/4)(m-l)]j 


(C2)1'2  =  4 


=  wavelength  at  spectrum  peak 


=  gamma  function 
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Spectral  Spreading  Functions 


Wrapped  Normal  Spreading  Function.  A  convenient,  easily  interpreted,  and 
mathematically  tractable  spreading  formula  is  provided  by  the  wrapped  normal  directional  density 
(Mardia,  1972), 


D(0;f) 


—  +  E 


cos[n(0  -0O)] 


-  E 


exp{-(0  -0O  -  2nqr/2cz} 
\[2n  c 


where  c  is  the  circular  standard  deviation  in  radians.  The  first  formula  is  best  if  c  >  n,  while 
the  second  is  better  if  c  <  n. 

For  the  usual  case,  where  c  <  <  ir/2,  the  half-peak  width  of  the  spreading  function  is: 

0HP  =  c/81^2  C64) 


Von  Mises  Function  (Mardia,  1972).  Other  common  models  from  directional  statistics 
are  the  von  Mises  formula  (here  10  =  modified  Bessel  function  of  order  zero): 


expfacos(0  -  0n)] 

D(9;f) - „  ,  * 

2  ic  ^(a) 


with  half-peak  width  of 


0HP  =  2  cos 


I" 


Generalized  Cosine-Squared  Function. 


0  _  0 

D(0;f)  -  K  cos2B  - - 

2 


with  half-peak  width 


0HP  =  4  cos'1  [(0.5),/2“] 
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A  reasonable  set  of  equivalent  values  for  c,  a,  and  a  are  obtained  by  equating  the  half-peak 
widths, 


_  4  cos'1  [(0.5)1/2*]  (69) 

and  then  solving  for  one  in  terms  of  the  other 
Water  Particle  Kinematic  Stretching 

Wave  properties  above  mean  water  level  are  not  really  defined  within  linear  wave  theory. 
Yet,  this  region  is  of  profound  importance  to  the  design  of  engineering  structures.  Various 
approximate  schemes  to  remedy  this  difficultly,  at  least  approximately,  have  been  suggested 
(Rodenbusch  and  Forristall,  1986). 

Functional  Extrapolation.  The  most  simple  procedure  would  be  to  use  the  formulas 
directly  from  linear  wave  theory,  with  z  given  values  above  mean  water  level.  If  z  is  measured 
positively  downward  from  mean  water  lever,  this  would  involve  inserting  negative  z  values  into 
the  formulas  for  the  kinematics.  The  derivation  for  linear  wave  theory,  however,  assumes  that 
0  <  z  <  d.  Asa  consequence,  this  approach  leads  to  predictions  of  velocities  at  the  crest  that 
are  much  too  large. 

Truncation.  Another  simple  procedure  is  to  use,  for  the  wave  kinematics  above  mean 
water  level,  the  value  of  those  same  kinematics  at  z  =  0.  However,  that  doesn’t  seem  to  lead 
to  reasonable  values  either. 

Linear  Extrapolation.  Linear  extrapolation  is  based  on  using  the  rate  of  change  of  the 
wave  property  with  respect  to  z,  at  z  =  0  and  linearly  extrapolating  for  values  above  mean  water 
level. 

Let 


c/8log2  «  2cos_1 


p0  =  wave  property  value  at  z  =  0 
p  =  dp(z)/dz  at  z  =  0 


(70) 


Then  the  linear  extrapolation  formula  is: 


p(z)  =  p0  +  p£z 


(71) 
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Reid-Wheeler  Stretching.  Another  scheme  that  has  had  some  use  in  published 
engineering  studies  (Wheeler,  1969)  consists  of  proportionally  stretching  the  wave  property  at 
z  =  0  up  to  z  =  -  T) ,  where  rj  is  the  water  level  elevation  above  mean  water  level  (i.e.,  positive 
upward).  Then,  the  linear  wave  theory  is  computed  with  the  vertical  coordinate  equal  to  zs, 
where: 


(d-z#)/d  =  (d  -  z)/(d  +  x\)  (72) 

Because  -ti  <  z  <  d,  it  follows  that 

zjd  =  1  -  (d  -  z)/(d  +  q)  (73) 


is  always  between  0  and  1 . 

Delta  Stretching.  Studies  by  Rodenbusch  and  Forristall  (1986)  suggest  the  following 
empirical  procedure,  which  seems  to  produce  fair  agreement  with  field  test  measurements.  Let 
0  ^  A  s  1.0  and  0  s;  dA  £  d  be  calibration  numbers  and  define  zA  so  that, 

(dA  -  zA)/(t|  A  +dA)  =  (dA  -  z)/(dA  +  n)  (74) 


If  a  =  0  and  dA  =  d,  then  the  Reid-Wheeler  stretching  listed  above  is  produced.  If  A  =  1  and 
dA  =  d,  then  zA  =  z  and  the  use  of  zA  would  correspond  to  functional  extrapolation. 
Rodenbusch  and  Forristall  suggest  the  following  empirical  rules: 

1.  If  dA  <;  z  <;  d,  use  z  directly  in  the  linear  wave  theory. 

2.  If  z  <  dA,  compute  zA:  (a)  if  0  £  zA  s  dA,  use  zA  to  compute  the  wave  kinematics: 
(b)  if  zA  <  0,  compute  the  wave  property  by  the  linear  extrapolation  as: 

p(z)  -  P0  +  P^zA  (75) 

After  studying  various  field  data,  Rodenbusch  and  Forristall  found  fairly  satisfactory  values  for 
A  and  dA  to  be: 


A  =  0.3 , 


06) 


where  an  is  the  standard  deviation  of  the  sea  surface. 
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Gamma  Extrapolation.  Another  formula,  suggested  here  as  purely  conjectural,  may  be 
based  on  the  gamma  function.  Suppose  the  kinematic  property  above  mean  water  level  is 
modeled  by: 


p,(z)  -  ce-(***wP(z  +  n)*'1 


(77) 


for  -n  <  z  <  0.  Let  p0  and  po  be  as  defined  for  linear  extrapolation.  If  the  gamma  curve  is 

matched  to  have  ps(z)  =  p  and  p,  =  Po ,  and  if  the  mode  of  the  gamma  is  forced  to  occur  at  z 
=  yrj,  where  0  <  y  <  1  ,  it  is  easy  to  define  a,  /3,  and  p(z): 


P 


(Y  ~  1)  Pq 

/ 

Po 


(78) 


p(z) 


exp(-z/  p) 


(79) 


a 


1  + 


n  Y 

P 


(80) 


This  curve  has  several  appealing  features  in  that  it  can  be  forced  to  attain  the  maximum  just 
below  the  water  surface,  but  is  zero  at  the  air-water  interface  (where  it  presumably  would  drop 
abruptly  to  zero).  Also,  it  matches  the  kinematics  at  z  =  0  and  has  the  matching  slope  there. 
However,  it  has  not  been  matched  against  real  data,  so  it  is  suggested  here  as  only  an  interesting 
possibility. 

Comparisons  with  Data.  Comparisons  of  the  delta  stretch  procedure  for  three  different 
data  sets,  one  from  a  wave  tank  study  and  two  from  field  measurements,  are  reported  by 
Rodenbusch  and  Forristall  (1986).  They  showed  fairly  good  empirical  correspondence,  although 
there  were  discrepancies.  Stretching  is,  at  best,  an  empirical  ad  hoc  procedure  and  should  not 
be  expected  to  yield  perfect  prediction. 


RANDOM  SIMULATION  OF  IRREGULAR  WAVES 

The  most  common  assumption  used  for  random  simulations  of  irregular  waves  is  that  the 
wave  properties  form  a  Gaussian  stochastic  process.  This  means  that  all  the  wave  properties  at 
any  selected  set  of  times  and  locations  follow  a  multivariate  normal  probability  law.  It  is  also 
usually  assumed  that  separate  components  are  independent  of  each  other  and  that  the  real  and 
imaginary  parts  of  the  complex  amplitudes  are  independent  of  each  other.  As  shown  in  the 
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Appendix,  the  real  and  imaginary  parts  of  the  amplitude,  Amj,  are  independent  and  normally 
distributed  if,  and  only  if,  the  0mj  are  uniform  on  (0,  2*),  the  a,^  are  Rayleigh  distributed,  and 
0mj  and  a„,j  are  independent  of  each  other. 

The  Rayleigh  random  variable  here  has  a  distribution  function: 


F(^p  =  I  ~  <*p(~  *l,/2  °L,) 


(81) 


where 


<i  -  S(fJD(0j;fJAfAd/2 


(82) 


A  convenient  way  to  generate  a  random  number  with  a  specified  probability  is  to  produce 
a  uniform  random  number  on  the  interval  (0,1)  and  then  to  equate  this  to  the  distribution  function 
for  the  desired  special  random  number  (Zelen  and  Severe,  1964,  p.  950).  Let  Uj  be  a  uniform 
random  number.  Then, 


F(amj)  =  U 


(83) 


gives 


amj  =  ^mj  v/-2  logC1  -  u,) 


(84) 


These  can  be  combined  to  obtain  a  random  simulation  of: 


Ki  =  amie 


(85) 


An  alternative  and  equally  good  approach  is  to  directly  produce  two  independent  standard  normal 
random  numbers,  Zj  and  Z Then, 


Ki  =  °*}<Zl+iZ2) 


(86) 


Either  method  can  be  used  to  build  up  a  simulation  of  the  complex  amplitude  matrix,  from  which 
all  the  linear  wave  properties  can  be  produced  through  the  application  of  the  appropriate  transfer 
functions. 

It  is  often  convenient  to  skip  the  Rayleigh  part  of  the  simulation,  set  a,^  =  am:,  and  use 
only  the  random  phase  part.  If  the  directional  spreading  is  not  too  narrow,  so  that  at  least  four 
or  five  energetic  directional  bands  are  added  at  each  frequency,  a  central  limit  theorem  behavior 
makes  the  resulting  wave  properties  behave  as  multivariate  normal  random  variables  with  the 
correct  spectra  and  cross  spectra;  this  is  called  the  random  phases  model  (Borgman,  1982b,  p. 
412).  However,  if  the  Rayleigh  part  is  retained,  a  11  components  will  exactly  introduce  Gaussian 
behavior.  An  example  of  wave  profile  simulation  is  given  in  Figure  3. 
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Figure  3 

An  unconditional  simulation  of  a  wave  profile  (variance  =  50,  k  =  1 .0,  fQ  =  0.1, 
principal  direction  =  60°  (from),  spreading  standard  deviation  =  20°). 


The  user  of  multidirectional  random  phase  wave  simulations  or  the  Rayleigh  amplitude 
wave  simulations  should  carefully  distinguish  between  the  input  spectral  density  (i.e.,  the 
theoretical  spectral  density  of  the  population)  and  the  sample  spectral  density  that  can  be 
computed  from  the  simulations.  The  sample  spectral  density  will  contain  ordinary  sampling 
variations  from  the  input  spectra,  and,  thus,  may  depart  rather  substantially  from  the  input 
function,  both  in  overall  variance  and  even  in  functional  shape.  If  the  sample  spectra  are 
computed  by  FFT  procedures,  the  deviations  between  the  sample  and  input  spectral  densities  will 
be  related  to  a  chi-squared  probability  density. 

Occasionally,  users  choose  to  simulate  with  the  random  phases  model  for  unidirectional 
wave  trains  because  this  has  the  advantage  that  every  simulation  has  its  sample  spectral  density 
exactly  equal  to  the  input  spectral  density.  Comparisons  of  platform  response  to  wave  forces  for 
such  simulations  then  show  the  variations  possible  between  seas  that  exactly  have  the  specified 
spectral  density.  These  simulations  are  examples  of  what  may  be  called  constrained  simulations. 
That  is,  the  simulation  is  forced  to  have  exactly  some  particular  behavior  (specified  spectrum  in 
the  unidirectional  random  phases  simulation  case)  that  it  would  not  have  in  a  freely  varying 
simulation. 

The  unidirectional  random  phases  simulation  has  the  additional  property  that  it  is 
approximately  Gaussian  in  the  time  domain  (unless  the  spectrum  is  very  narrow-banded)  through 
a  central  limit  theorem.  However,  it  is  not  Gaussian  in  the  frequency  domain,  because  the  FFT 
coefficients  are  of  the  form  Am  =  c  exp(-it),  where  c  is  a  deterministic  constant  and  *  is  a 
uniform  random  variable  on  (0,  2rr). 

A  multidirectional  random  phases  simulation  is  approximately  Gaussian,  even  in 
frequency,  because  now  Am  =  2  Cj  exp(i*j)  and  the  summing  over  directions  is  enough  to  cause 
approximate  normality,  at  least  in  the  nonextreme  ranges  of  simulated  data,  through  a  central 
limit  theorem.  A  central  limit  theorem  formulation  suitable  to  these  derivations  is  given  by 
Takano  (1954)  in  a  doubly  subscripted,  multivariate  format.  A  typical  derivation  of  the  Gaussian 
properties  of  random  phase  models  is  provided  by  Brown  (1967). 

Another  property  of  simulations  that  should  be  understood  by  users  is  that  of  ergodicity. 
In  words,  a  process  that  is  ergodic  should  have  its  probabilistic  behavior  for  infinitely  long  time 
intervals  be  the  same  as  the  probability  laws  for  the  ensemble  of  possibilities  at  a  single  time  or 
at  several  specified  times.  Because  the  simulation  procedures  presented  previously  are  based  on 
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a  fixed  matrix  of  complex  amplitudes,  Amjl  determined  once  and  for  all  by  a  finite  set  of  random 
numbers,  the  time  series  given  by  the  Amj  will  always  have  the  random  bias  introduced  by  the 
random  number  sampling  variability,  regardless  of  how  long  the  time  series  is  made.  Thus,  the 
simulation  process  provided  is  not  ergodic  as  N  tends  to  infinity  (see  Jefferys,  1987;  Tucker,  et 
al.,  1984). 

Tucker  is  particularly  concerned  that  the  standard  deviation  of  the  wave  property  be  the 
same  from  (x,  y)  location  to  (x,  y)  location,  because  such  behavior  is  important  in  his 
applications.  In  real  wave  data,  if  the  time  interval  of  sampling  is  long  enough,  and  if  the  wave 
field  statistics  remain  stationary  over  the  multiple  hours  of  record  required,  ergodicity  will  force 
the  variance  at  every  spatial  location  to  be  the  same.  However,  in  the  usual  20  minutes  of  most 
wave  data,  there  will  be  substantial  differences  in  variance  from  location  to  location.  Real  world 
data  show  the  same  spatially  varying  sampling  variability  as  the  simulations. 

Simulations  of  directional  seas  on  the  multidirectional  case  can  be  constrained  to  have 
exactly  a  specified  frequency  spectral  density.  The  procedure  is  as  follows.  The  matrix  of 
complex  amplitudes  Amj  is  simulated  with  full  sampling  variability  as  shown  in  Equation  86. 
Then  the  sequence  (Amj;  j  =  1,...,J}  is  summed  as: 

E  IA.il2 

n  =  ii! _  (87) 

E  < 

j.i 


Finally,  the  original  Amj  are  divided  by  B1/2  to  get  a  new  set  of  Amj  at  that  frequency.  This  is 
repeated  at  the  other  frequencies  so  that  the  simulation  will  have  exactly  the  specified  S(mAf) 
spectral  density. 

A  similar  procedure  can  be  used  to  constrain  the  directional  simulation  to  have  a 
directional  spectral  density  that  exactly  agrees  with  the  theoretical  population  directional  spectral 
density.  Here  the  complex  amplitude  matrix  is  simulated  as  before  with  complete  sampling 
variation  freedom.  Then  the  Amj  =  Umj  -  iVmj  are  replaced  by: 

A»j  %/(<  *  <t)m  (88) 


The  important  thing  to  remember  is  that  a  simulation  is  always  artificial.  It  will  retain  some 
properties  of  the  real  world  but  will  violate  others.  There  is  not  one  correct  way  to  compute 
simulations,  but  various  ways  that  are  appropriate  to  particular  applications.  Jefferys  (1987)  and 
Tucker,  et  al.  (1984)  provide  a  real  service  in  emphasizing  the  dangers  in  blindly  using  a 
simulation  procedure  without  appreciating  the  artificial  aspects  of  the  simulation.  However,  one 
should  not  go  too  far  the  other  way  and  condemn  a  simulation  procedure  for  all  applications, 
when  it  is  perfectly  satisfactory  in  many  problems  where  those  particular  defects  are  unimportant. 

Usually,  ergodicity  is  not  particularly  important  for  studies  of  wave  forces  on  offshore 
structures  as  long  as  the  sampling  variability  for  a  given  finite  interval  of  simulation  is 
appreciated.  A  typical  treatment  involves  looking  at  the  force  behavior  in  a  number  of  separate 
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simulations.  This  represents  a  sample  "ensemble  set"  and  convergence  to  the  correct  probability 
laws  or  correct  spectral  density  will  occur  by  averaging  over  an  increasing  number  of  simulation 
sets. 

A  very  interesting  application  of  the  unconditional  simulation  procedures  given  in  the 
preceding  pages  has  been  developed  by  M.  J.  Briggs  at  the  Coastal  Engineering  Research  Center, 
Vicksburg,  Mississippi  (Briggs,  et  al.,  1987).  He  simulates  water  level  elevation  time  histor^ 
at  directional  paddle  locations  with  an  input  directional  spectral  density,  then  converts  these  time 
series  to  stroke  magnitude,  and  finally  produces  a  directional  wave  field  within  the  wave  basin. 
He  then  reverses  the  procedure  by  using  an  array  of  model  tank  wave  staffs  in  the  center  of  the 
basin  as  input  to  a  directional  spectrum  analysis  program  and  estimates  the  directional  spectral 
density  that  has  been  developed.  Interestingly  enough,  the  estimated  spectra  agree  quite  well 
with  the  original  input  spectra,  within  the  usual  range  of  sampling  and  estimate  variability. 


CONDITIONAL  SIMULATION  OF  IRREGULAR  WAVES 

The  various  linear  properties  of  the  irregular  wave  train  will  be  assumed  to  follow  the 
multivariate  normal  probability  law  over  space  and  time.  Within  this  framework,  conditional 
probability  laws  are  themselves  multivariate  normal  and  an  elaborate  theory  can  be  constructed. 
It  will  be  assumed  in  the  following  that  the  reader  1m  some  familiarity  with  multivariate 
statistical  analysis  and  matrix  theory. 

The  problem  that  will  be  considered  here  is  the  following.  Suppose  wave  properties  have 
been  measured  for  some  set  of  wave  conditions  of  interest,  and  the  directional  spectrum  S(f,0) 
has  also  been  estimated  or  an  acceptable  model  has  been  found.  What  are  reasonable  time  series 
for  other  wave  properties  co-occurring  with  the  measured  set,  or  for  the  same  wave  properties 
at  times  other  than  the  measurement  intervals? 

The  wave  properties  that  were  not  measured  are  stochastic  processes  that  are  to  some 
degree  correlated  with  measured  data.  Depending  on  the  extent  and  amount  of  the  actual  data, 
the  nonmeasured  wave  properties  may  be  constrained  to  strong  agreement  with  the  measurements 
or  may  only  be  weakly  related  to  them. 

Conditional  simulation  of  wave  property  time  series  statistically  consistent  with  a  specified 
measurement  set  provides  a  very  powerful  approach  to  certain  ocean  engineering  problems.  The 
method  does  not  provide  a  single  determination  of  the  nonmeasured  time  series,  but  rather  one 
for  each  simulation.  The  range  of  variation  from  simulation  to  simulation  provides  a  measure 
of  how  strongly  the  measurements  determine  the  structure  behavior. 

Although  these  methods  are  quite  new  in  ocean  engineering,  they  have  had  some  use  in 
the  petroleum  industry  (Rodenbusch  and  Forristall,  1986)  and  the  author  is  aware  of  at  least  two 
other  studies  where  these  methods  were  highly  successful  in  applications  to  ocean  waves 
(Vartdal,  Krogstad,  and  Barstow,  1989).  Rather  similar  methods  in  geostatistics  and  mining 
(Joumel,  1974;  Borgman,  Taheri,  and  Hagan,  1983)  have  been  used  for  spatial  random  fields 
to  some  extent  for  several  decades. 

The  usual  computer  simulation  of  waves  (Borgman,  1982b)  satisfying  a  specified  model 
for  the  directional  spectral  density  suffers  from  a  serious  practical  defect  if  one  is  primarily 
interested  in  producing  very  large  waves.  Most  simulations  produce  average  waves.  Very  long 
computer  runs  are  required  to  capture  the  occurrence  of  an  extra  large  wave.  However,  the 
conditional  simulation  starts  with  the  inclusion  of  a  large  wave  profile  embedded  into  the  wave 
train.  The  various  associated  kinematics  are  produced  consistent  with  this  large  wave,  for  each 
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simulation.  Every  computer  run  produces  another  simulation  that  is  appropriate  to  the 
engineering  analysis. 

A  large  variety  of  problems  can  be  approached  from  the  conditional  simulation 
perspective.  In  addition  to  using  measured  storm  wave  intervals  to  condition  on,  one  could 
introduce  wave  groups  of  three  or  four  large  waves  in  sequences,  condition  on  these  to  produce 
the  associated  kinematics,  and  then  impose  all  this  onto  a  structure  for  dynamic  behavior. 
Another  type  of  investigation  might  be  to  study  the  variability  in  kinematics  associated  with  a 
large  Stokes’-type  wave  profile,  when  a  directional  spectrum  is  present  with  various  degrees  of 
spreading.  Undoubtedly  many  more  diverse  applications  will  be  found  as  these  methods  become 
more  widely  known. 

Multivariate  Normal  Probability  Law 

The  n  component  vector  V  with  mean  p  and  covariance  matrix  C,  will  be  said  to  behave 
according  to  a  multivariate  normal  probability  law  if  its  probabilities  of  occurrence  follow  the 
probability  density: 


Fv(v)  =  [(2  it)*12  /jcT ]  1  exp[-(v  -  p)TC_I  (v  -  p)/2]  (»9) 

where  the  superscript  T  denotes  the  transpose,  |C|  is  the  determinant  of  C,  and  C'1  is  the 
inverse  of  C.  The  inverse  of  C  in  the  exponent  can  be  replaced  with  a  generalized  inverse,  C  +  , 
if  suitable  restrictions  to  the  domain  space  are  introduced  (Pringle  and  Rayner,  1971 ,  pp.  70-72). 
It  is  often  more  convenient  to  use  a  shorthand  notation  for  the  specification  of  multivariate 
normality  as: 


V:  N(p,C) 


(90) 


This  is  read  as  V  is  multivariate  normal  with  mean  p  and  covariance  matrix  C. 

Conditional  Normal  Density 


Suppose  the  random  vector  V  can  be  broken  up,  or  partitioned,  into  two  parts,  Vj  and 
V2,  and  that  the  combined  vector  satisfies: 


V 


fv, 


:  N 


Pi 

P2 


c  c 

^11  v'12 

c  c 

'-21  v-22 


(91) 


That  is, 
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E[Vjl  -  |i, 

E[V2]  -  p2 

covariance  matrix  of  (Vj)  =  E[(Vj  -  p,)(Vj  -  p,)T]  =  Cn 
covariance  matrix  of  (V2)  =  Cn 

cross-covariance  matrix  of  (Vj,V2)  =  E[(Vt  -  p,)(V2  -  p2)T]  -  C12 


(92) 

(93) 

(94) 

(95) 

(96) 


Suppose  specific  values  for  Vj  are  assigned  numerically  as  Vj  =  Vj.  Here  Vj  is  a  list  of 
numbers,  while  Vi  is  the  abstract  random  vector.  Then  the  probability  density  for  V2,  given  Vj 
=  Vj,  is  (Anderson,  1958): 


(V,|V,  *  v,):  N(i4  ♦cjq-,,(v1  -  It,).  CB-c5cr,'CB) 


(97) 


In  its  simplest  form,  conditional  simulation  is  just  any  procedure  for  producing  a  vector  V2  with 
the  mean  vector  and  covariance  matrix  given  above. 

Conditional  Simulation  Method 

A  straightforward  method  (Borg man,  Taheri,  and  Hagan,  1983)  is  given  by  a  two-step 
procedure: 

1.  Unconditionally  simulate  V,  without  reference  to  Vj  =  v,.  Let  Vj  and  V2  denote 
these  unconditional  simulations. 


2.  A  conditional  simulation  is  provided  by: 


(VjIV.-v.)  -  cJc,V(»1-v1)+vJ 


(98) 


All  of  the  techniques  based  on  this  procedure  require  substantial  facility  in  deriving  and 
computing  the  covariances  in  Cn  and  C12  and  in  inverting  Cu.  For  example,  consider  two 
general  wave  properties  of  the  form  of  Equation  25  combined  with  Equation  39.  Let  X  —  1,2 
[  with  At  =  time  increment,  Af  =  l/(NAt)]: 


B2(m  Af;x,y,z)  =  A  =  G2(m  Af,z)  T2(m  Af)  HA(j  A 6) 

H  1  (99) 

*  exp [-i  p0 l^fx  cos(j  A0  -  0X)  +  y  cos(j  A0  -  0y)J  ] 

for  0  <  m  <  N/2.  Set 
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Ba(N  Af/2;x,y,z)  =  0 


(100) 


BX([N  -m]  Af;x,y,z)  »  Bx  (m  Af;x,y,z) 
for  N/2  <  m  <  N.  The  wave  properties,  in  time  domain,  are  given  by: 

N-l 

bx(n  At;x,y,z)  =  £  Bx(m  Af;x,y,z)exp(-i2  n  mn/N)  l]01) 

m-0 


for  n  =  0,1,2,...,N-1. 

The  covariances  between  two  time  domain  wave  properties  are  developed  in  Theorem  D 
of  the  Appendix.  That  is: 

QiXl  =  Gx(mAf,z)Tx<mAf)Hx(jA0) 

*  expj-ipokjx  cos(j  A0  -  0X)  +  y  cos(j  A0  -0y)J  ]  (102) 

Equation  2  of  the  theorem  gives  the  covariance  between  the  two  wave  properties  one  at  time  n 
At  and  the  other  at  time  n'  At.  The  location  coordinates  (X|,yj,Zj)  for  the  first  wave  property 
must  be  substituted  into  OS’  while  the  corresponding  values  for  location  of  the  second  wave 

property  (x^y^z^  must  be  substituted  in  Q®.  Stated  in  terms  of  integral, 


-  2  x 

E[b(1)(t)b®(t  +  t)]  =  /  /  S  (f, 0)  Q(1)*  (f, 0)  Q®(f, 0)  exp(i  2  x  f  x)  d0  df  (103) 

—  o 

The  computation  of  the  wave  property  covariances  can  be  quite  a  task.  Some 
simplification  and  acceleration  of  the  calculations  can  be  achieved  by  introducing  the  Fast  Fourier 
Transform  algorithm.  However,  this  still  leaves  a  substantial  effort  in  manipulating  Cj2  and 

computing  cn'  for  the  often  large  matrices  that  occur  in  realistic  applications. 

Conditional  Simulation  of  Complex  Wave  Amplitudes 

Great  savings  in  computer  run  time,  as  well  as  the  achievement  of  very  substantial 
simplification  and  clarification  of  the  model  in  application,  is  obtained  by  centering  attention  on 
the  matrix  of  complex  wave  amplitude,  Amj  =  Umj  -  iVmj  introduced  previously.  The  Amj  are 
conditionally  simulated,  given  the  specified  wave  properties  and  the  selected  directional  spectral 
density. 


Figure  4 

Two  conditional  simulations  of  water  particle  velocities,  given  the  profile  in 
Figure  3.  (Waves  from  heading  of  60°:  x  axis  east:  y  axis  north: 
velocities  at  15.7  m  below  mean  water  level.) 


The  actual  calculation  of  these  complex  amplitude  simulations  requires  a  fairly 
sophisticated  introduction  of  frequency  domain  methods  as  interrelated  to  time  domain  data.  The 
Appendix  provides  a  very  succinct  summary  of  a  fundamental  set  of  theorems  giving  a  theoretical 
framework  for  the  discrete  Fourier  transform. 

Frequency  Domain  Conditioning 

A  common  special  case  in  applications  arises  when  the  wave  profile  or  other  wave 
properties  have  been  measured  over  the  full  time  interval  of  interest.  Perhaps  this  represents  an 
interval  of  time  during  a  historical  storm  when  wave  action  was  quite  high.  A  simulation  of  the 
wave  kinematics  throughout  the  wave,  conditioned  to  agree  with  the  occurrence  of  the  measured 
time  series,  is  needed  for  wave  force  analysis.  What  is  the  most  efficient  way  to  develop  this 
conditional  simulation? 

A  very  effective  technique  is  to  first  perform  FFT  on  the  measured  time  series  to  obtain 
the  Fourier  coefficients,  and  then  to  conditionally  simulate  each  Amj  =  Umj  -  iVmj  separately, 
given  the  Fourier  coefficients  {Bx(m  At;  x,  y,  z),  X  =  1,...,  number  of  measured  time  series} 
for  the  measured  wave.  This  will  be  illustrated  with  the  case  for  one  measured  time  series  and 
the  case  for  two  measured  time  series.  Figure  4  gives  examples  of  such  simulations. 
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Consider  first  the  situation  for  one  measured  time  series  {b(n  At),  n  =  0,1,..., N-l},  with 
FFT  coefficients: 


N-l 

B(mAf)  =  £  b(n  At)exp(-i2  n  mn/N)  00*) 

B-0 

with  the  Q,^  from  Equation  102,  B(mAf)  can  be  related  to  the  complex-valued  wave  amplitude 
as: 


B(m  Af)  =  £  A*jQaiJ 


i-i 


(105) 


The  goal,  then,  is  to  conditionally  simulate  Umj  and  Vmj,  give  B(m  Af)- 

The  following  structural  assumptions  will  be  made: 

1.  All  time  series  will  be  assumed  to  be  covariance-stationary  with  zero  mean. 

2.  All  time  series  will  be  repeated  periodically. 

3.  Aqj  =  An/2ij  =  0  for  all  j  =  1,2,...,J. 

The  imposition  of  the  second  assumption  has  the  effect  of  making  the  last  part  of  the  simulation 
have  a  correlation  with  the  first  part.  If  this  causes  improper  behavior  in  the  application  of  the 
series,  about  three  wave  periods  of  simulation  should  be  deleted  at  the  end  of  the  simulated  time 
series,  before  it  is  used.  This  just  means  that  N  is  made  a  little  larger  than  is  needed  in  the 
application  to  adjust  for  the  deletion  of  some  time  steps  at  the  end  of  the  simulation.  The 
addition  of  the  periodicity  assumption  produces  such  a  clean,  elegant  statistical  structure  in  the 
Fourier  coefficients  that  it  is  definitely  worth  the  trouble  of  compensating  for  the  minor 
distortions  introduced  by  the  assumption.  The  imposition  of  the  third  assumption  just  forces  the 
simulations  to  have  mean  zero  and  to  have  no  energy  at  the  Nyquist  frequency. 

The  covariance  matrix  for  Amj  and 

B(m  Af)  =  +  i  Ym  (106) 


can  be  constructed  from  Theorem  D  in  the  Appendix.  It  is  best  stated  in  terms  of  the  real 
random  variables  Umj,  Vmj,  *m  and  Tm.  For  0  <  m  <  N/2: 
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Cov 
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*. 

U.J 


m 


ES.jIQ.jI2  0  S.^Q.,)  S^ImCCLp 

J-i 


0  E  S.j|Q.jf  S.jInKQ.p  -S^RctQ.p 

S»jRe[Qmji  SMjlm(Qmj)  S.,  0 

Smj^m(Q.j)  -s^Reio^]  0  sfcj 


007) 


The  actual  conditional  simulations  are  then  just  an  application  of  Equation  98,  with: 


C 


u 


E  S^IQ^fAfAe^ 

J-l 


(108) 


c 


12 


Re(Qmj)  toCQmj)  Smj  Af  A0 
^(Qmj)  -^((^pj  2 


(109) 


V 


l 


from  the  measured  time  series 


(HO) 


The  unconditional  simulation  of: 


(u  • 

mj,ui 

v 


(HD 


can  be  produced  easily  from  two  independent  standard  normal  random  numbers  (Zj  ,Z^)  with  the 
formulas: 


Uaj>».  -  Z^AfAe/l)1* 


(112) 
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ZjCS.jAfA  0/2)1/J 


(113) 


This  is  repeated  for  the  other  j-indices  with  additional  pairs  of  j-indices  with  additional  pairs  of 
independent  standard  normal  random  numbers.  The  unconditional  simulations  of  *m  and  i|rm  &e 
completely  determined  by  the  set  of  values  {(U^  ^  j  =  1,2,...,J}  as: 


Q  *  iY  =  y  (U  .  -iV  .  )0 

m,ui  m,ui  Z-f  '  «nj,u»  mj,u*  '  V£»j 

j*l 


(114) 


Equation  98  for  the  conditional  simulation  becomes: 


(Val^-v,)  = 


V-W 


[Re(Qmj)  Im(Qmj) 


SmiAfA0/2 


[lm(Qmj)  -Re(Q.pJ  E  SmjjQrajf  AfA0/2 


«  1  f«  1 
“  ».»* , 


(115) 


The  validity  of  Equation  1 15  can  be  verified  by  confirming  that: 


®m.cs  ~  E  ^mj.cs  "  1  \nj,c*)  Qmj 
j-1 


(116) 


is  the  same  as  B(m  AO  in  Equation  106.  That  is,  the  complex  amplitudes  produced  by 
conditional  simulation  do,  in  fact,  produce  the  FFT  coefficients  for  the  conditioning  time  series. 
The  proof  of  Equation  116  follows  from: 
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(117) 


d,-i)  t-ic- 

Vv»j.c») 


Applying  this  to  Equation  1 15,  one  gets,  after  some  algebra, 


Q^j  SB^Af  A0 


i:  s.j|Q.,f4f4e 

j-1 


[B--B«.«]  +  A.].. 


(118) 


The  substitution  of  Equation  118  into  Equation  116  gives: 


B  =  B  -  B  +  V  A.  Q_.  =  B 

m.ci  m  m,u*  Taj.M^mJ  n 

j-1 


(119) 


which  is  the  verifying  equality. 

The  covariance  is  more  complicated  if  there  are  two  time  series  that  are  simultaneously 
measured,  which  are  to  be  used  for  conditioning.  Then, 

B?  =  E  Kt<&]  -  (120> 

j-i 


b?  -  E  Ki<£\  -  •?+**?  <121> 

j-i 


are  the  conditioning  FFT  coefficients.  The  covariance  matrix  which  must  be  evaluated  is  that 
for  the  six  variables  v®  a®  Yp)  U  ,.V  ,)•  For  this, 

i  A m  > *  *m  >  mj ’  mj' 
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The  conditional  simulation  proceeds  in  two  steps  as  before.  First,  an  unconditional  simulation 
of  Umj  and  Vmj  is  formed  with  pairs  of  independent  standard  normal  random  numbers  {Zmjl, 
Zmj2;  0  <  m  <  N/2,  j  =  with  the  formulas: 

IV.  -  Zmll  (Sml  Af40/2)1«  (124) 


)»  025) 


The  unconditional  simulation  of  and 


is  achieved  with  Equation  105  as: 


B®.  =  E(U„,.„-iV„J,u.)Qi1’ 


0) 


j-1 


(126) 
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(127) 


E  <«.!.« -*v.,.„)QS 

J-l 


The  second  step  consists  of  substituting  all  these  results  into  Equation  98.  It  can  be  shown,  as 
before,  that  the  conditional  simulation  of  the  A^,  when  used  with  the  transfer  functions  and 

q2},  reproduce  the  FFT  coefficients  used  in  the  conditioning.  However,  the  details  are  too 
complicated  for  presentation  here. 

Time  Domain  Conditioning 

Time  domain  conditioning  is  addressed  in  the  following  problem.  Suppose  one  or  more 
time  series  are  measured  over  portions  of  the  time  interval  of  interest.  The  problem  is  to 
simulate  conditionally  the  time  series  for  the  rest  of  the  time  interval  and  for  other  wave 
properties  during  the  entire  interval.  Figure  5  provides  an  example  of  this  technique. 


Figure  5 

Two  conditional  simulations  of  the  wave  profile  for  30  <  t  <  55, 
given  the  profile  in  Figure  3  for  0  <  t  <  30. 
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The  technique  will  be  illustrated  with  a  single  conditioning  time  series,  consisting  of  the 
water  level  elevations,  n( n  At).  Suppose  {r?(n  At),  n  =  0,l,2,...,v},  where  v  <  N  is  known 
and  is  to  be  used  for  conditioning.  The  covariance  matrix  for  {n(0),rj(2  At),»7(2At),...,i?(N0 
At),Umj,Vmj}  forms  the  basis  for  the  simulation. 

Let 


N-J  J 


ck  =  E  E  Smjexp[i2  7ikm/N]A0Af 


(128) 


m-0  j»l 


Then  Cn  is: 


co  ci 


ci  co  ci 


^ll  -  C2  ci 


(129) 


Cv  Cv-1  Cv-2 


and  the  C12  matrix  is: 


S  .  Af A0 

mj 


C,2  = 


Smjcos(2  ii  m/N)  Af  A0 
Smj  cos(4  n  m/N)  Af  A0 
Smj  cos(6  it  m/N)  Af  A0 


Smjsin(2nm/N)  Af  A0 
Smj  sin(4  n  m/N)  Af  A0 
S  jStn(6  n  m/N)  Af  A0 

in  j 


(130) 


Smjcos(2  it  mv/N)  Af  A0  Smjsin(2  itmv/N)  Af  A0 


The  conditional  simulation  proceeds  by  first  using  Equations  124  and  125  to  get  Umj  us  and 
Vmj  us.  Then  the  unconditional  simulation  of  the  rj(n  At)  is  computed  from: 


N-l  I 


*lus(nA0  =  j;  E  AffljuIexp(i2nmn/N) 

m*0  j«l 


(131) 


Finally,  Equation  98  is  used  to  get  the  conditioned  set  of  complex  wave  amplitudes  with: 


n(0) 

bM(0) 

ri(Af) 

n„(AO 

V1  - 

1,(2  Al) 

.  vi  " 

nM(2  ao 

n(v  ao. 

i,M(v  AO 

(132) 


The  inversion  of  Cn  can  be  done  with  special  techniques  for  Toeplitz  matrices  (Golub 
and  Van  Loan,  1983,  pp.  125-133)  or  with  conjugate  gradient  iteration  (Golub  and  Van  Loan, 
1983,  pp.  362-369).  Any  technique  that  solves  Cnx  =  b  for  x,  given  b,  can  be  used.  If  Cn 

is  singular,  the  generalized  inverse  of  Cj  j  (Searle,  1983,  pp.  212-226)  can  be  used  in  place  ofCu' 

in  Equation  98. 


TECHNIQUES  FOR  COMPLIANT  STRUCTURES 

The  techniques  for  simulation  presented  previously  are  based  on  a  fixed  coordinate 
position  (x,  y,  z).  If  the  structure  is  moving,  there  are  difficulties  in  following  the  structure 
from  one  position  to  another.  One  way  to  proceed  would  be  to  compute  time  series  at  each 
intersection  for  an  extensive  three-dimensional  grid  over  the  (x,  y,  z)  space.  Then  the  velocities 
and  accelerations  at  any  arbitrary  space  location  could  be  interpolated.  However,  for  a  typical 
structure,  this  procedure  requires  such  formidable  amounts  of  computer  storage  as  to  be 
impossible  on  even  large  computers. 

An  alternative  procedure  is  presented  here  that  surmounts  this  difficulty  and  requires  only 
moderate  computer  capacity.  The  procedure  is  based  on  the  need  for  wave  kinematics  and  other 
wave  properties  only  within  a  relatively  narrow  rectangle  of  horizontal  coordinate  position  as 
compared  to  the  wavelengths  involved.  The  approach  involves  the  expansion  of  the  wave 
properties  in  Legendre  orthogonal  polynomials  in  x,  y,  and  z. 

Let 


P„00  =  b0  +  bj  x  + ...  +  bn  x  " 

Pn  (x)  =  b0*  +  b|*  X  + ...  +  b  *  X  * 


(133) 


be  the  Legendre  and  shifted  Legendre  orthogonal  polynomials  of  order  n  (Hochstrasser,  1964). 
The  coefficients  are  selected  so  that: 
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if  m  *  n 
if  m  —  n 


(135) 


The  Legendre  polynomials  are  defined  for  -1  <  x  <  1  and  the  shifted  Legendre  polynomials  are 
defined  for  0  <  x  <  1. 

The  first  several  Legendre  polynomials  are: 


P0(x)  -  1 

(136) 

Pi  00  =  X 

(137) 

p2(x)  -  (3x2  -  l)/2 

(138) 

p3(x)  -  (5  x3  -  3  x)/2 

(139) 

p4(x)  -  (35  x4  -  30x2  +  3)/8 

(140) 

p3(x)  -  (63  x5  -  70  x3  +  15  x)/8 

(141) 

The  shifted  Legendre  polynomials  are  related  to  the  regular  Legendre  polynomials  by  the 
relation: 


Pn’W  -  P„(2*-1) 


(142) 


The  use  of  orthogonal  polynomials  to  approximate  an  arbitrary  function  g(x)  defined  on 
(-1,1)  can  be  illustrated  as  follows.  Suppose  the  approximation  to  be  used  is: 
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g(x) 


N 


£ 


attPB(x) 


The  coefficients  a,j  are  chosen  from  a  "least-squares"  criterion.  Let  a,j  be  those  values  that 
minimize: 


Q  = 


/ 


N 


gto  -  E  a.pDw 

■  •0 


2 

dx 


(144) 


Then, 


aQ 


-  /2 


-l 


N 


g(x)  -  £  aapn(x) 
n  »0 


pk(x)dx 


(145) 


Q  is  at  an  extreme  if  dQ/d^  =  0  for  all  k  =  0,1,2, ...N.  This  reduces  to: 


n  1  i 

E  a*  /  PatoPkOO*1  =  /  g(x)PkOOdx 


n-0  ,| 


(146) 


But,  by  the  orthogonality  relation,  this  further  reduces  to: 


=  /g(x)pk(x)dx 


(147) 


The  similar  development  for  shifted  Legendre  polynomials  gives: 


ak‘  =  (2k +  1)  jf  g(x)pk'(x)dx 


(148) 


How  can  these  relations  be  applied  to  ocean  wave  kinematics?  The  essential  canonical  form  for 
a  linear  wave  property  is  given  in  Equation  23.  It  should  be  noted  that,  in  every  case,  G(z)  is 
either  1.0  (water  level  elevation)  or  is  of  the  form: 
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+  Sj  e 


1  +  Sj  e 


2k.d 


where  k,,,  is  the  wave  number.  Thus,  the  general  wave  property  b(n  At)  can  be  expressed  from 
Equation  35  as: 


N-l 

b(n  At)  =  £ 

m>0 


(150) 


where  Bm  is  the  FFT  coefficient  given  by: 


Bm  =  E  COS(0  -0X)  +y  COS(0  -0  )}  1 

j-i 


(151) 


for  sea  surface  elevations,  and  by: 


J  le'^'  +  s  e  k«(2d-*)j 

E  AmiTmBi  ~ - 1 - ~ 

*  exp[-i  P0km{x  cos(0  -  0X)  +  y  cos(0  -0y)}  ] 


(152) 


for  the  other  wave  properties.  These  can  be  expressed  as  a  sum  of  products  of  separate  functions 
of  x,  y,  and  z  as  follows: 

1 .  Sea  Surface 

j 

B»  •  E  Ki  exP  p  P0  K x  c0^0  -  o.) ) exP  p  P0  K y  '  ey)  ]  (153) 


2.  Other  Wave  Properties 
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Suppose  it  is  desired  to  obtain  a  good  representation  locally  in  the  vicinity  of  the  structure. 
Consider  the  volume, 


*0‘D1 

S  X 

s  Xo  +  D, 

(155) 

y0-D2 

*  y 

*  y0  +  D2 

056) 

0 

£  Z 

i  d 

057) 

(Here  z  has  been  taken  as  positive  downward  and  zero  at  mean  water  level.) 
It  is  natural  to  scale  function  as  follows: 

1 .  Sea  Surface 


Bm  =  E  COSC0  "e,)]cxP 
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(158) 


2.  Other  Wave  Properties 
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where  kg  is  a  selected  single  reference  wave  number.  For  many  applications,  the  second  term, 

e-k»(2d-*)  (163) 


is  negligible  because  depth  is  large.  For  the  moment,  suppose  that  this  second  term  can  be 
ignored  (it  will  be  reintroduced  later).  Then,  if 


u  =  (x-x^/ D, 

v  -  (y-y o)/d2 


(164) 
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we  have  the  following: 
1 .  Sea  Surface 


B«  -  E  cos(0-0,)+yocos(0-0,)]} 
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(165) 


2.  Other  Wave  Properties 
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If  these  are  substituted  into  Equation  150: 


N  N  N 

b(nAt)  =  E  E  E  Q,,p.y(nAt)p«(u)pp(v)pY*(w1) 

*•0  P  »0  t  »0 


(167) 


where  (sea  surface  case): 


N-l 


=  E 


m-0 


E  ®b  ^p  Am  j  ®*P  {  ~  *  P0  *0  COS(0-0,)} 


j-1 


*  exp{-i  P0kmy0  cos(0  -  0y)}]ei2’““'N 


(168) 


Equation  168  is  an  FFT  of  the  quantity  within  []  for  each  combination  of  a,  0,  and  y. 
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The  other  wave  properties  have  a  similar  expression,  with 


N  - 1 
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E‘.VA,T.H. 
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At  a  given  time,  nAt,  many  of  the  Qa>p  Y  for  a  given  wave  property  are  negligible. 
After  all,  the  region  x0  ±  Dj  and  y0  ±  D2  is’relatively  small  compared  to  the  wavelengths. 
Hence,  low  order  polynomials  are  all  that  are  required  in  order  to  represent  the  variation  over 
the  horizontal  region.  The  vertical  variation  is  attenuated  more  or  less  exponentially  with  depth, 
so  a  polynomial  in  w(  should  only  need  relatively  low  order. 

Hence,  at  a  given  time,  only  a  few  Qa,p)Y  will  be  needed  to  represent  the  wave  property. 
The  particular  coefficient  needed  may,  however,  be  different  from  one  time  step  to  another.  An 
interesting  point  here  is  that  the  Qa  «  Y  (nAt)  can  be.  themselves,  computed  by  the  FFT 
algorithm  simultaneously  for  n  =  0,1,2,..., N-l. 

Up  to  this  point,  the  actual  computation  of  aB ,  bp  and  cy  has  not  been  explicitly 

stated.  From  the  definition  of  orthogonal  polynomials, 

a.  =  JexpJ-ik^PoDj  cos(0  - ex)u]p.(u)du  (I70) 

2  -i 


BP  "  /e*P[-iMoD2C<*(0-0y)v]Pp(V>dv 
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cy  =  (2y  +1)  /w^p^WjJdw, 
o 


(172) 


Other  depth  term 
Let 
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(173) 


w2  -  e-*a--*) 


Then  an  exactly  similar  expansion  with  the  same  coefficients  can  be  developed.  The  resulting 
representation  of  the  wave  property  is: 

N  N  N 

b(nAt)  =  E  E  E  Q«pT(“AOp.(u)P,j(v) 

• *o  P-o  y-o  (174) 

|Py  (Wl>  +  Si  Py  (W2>  } 


Statistics  of  Random  Wave  Geometry 

In  testing  the  statistics  of  wave  properties,  such  as  height  and  period,  it  is  important  to 
distinguish  between  properties  for  individual  waves  and  average  properties  above  a  specified  time 
interval.  The  individual  wave  properties  are  the  properties  that  a  single  wave  possesses  at  a 
particular  instant  in  time.  Individual  wave  heights  are  usually  defined  in  terms  of  maxima  or 
minima  between  zero  crossings  of  the  water  level  elevation  above  mean  water  level  or  in  terms 
of  a  time-varying  envelope  function.  Individual  wave  periods  may  be  defined  in  terms  of  time 
interval  between  up-crossings  or  in  terms  of  functions  of  time  derivatives  of  various  orders  of 
the  water  level  elevation.  In  contrast,  time  interval  average  properties  such  as  significant  wave 
height  or  predominant  wave  period  are  usually  taken  from  the  wave  spectral  density  holding 
during  that  time  interval. 

There  is  a  very  extensive  literature  on  statistics  of  wave  properties.  Review  and 
expository  papers  are  given  by  Borgman  (1982a)  and  Ochi  (1982).  In  addition,  there  are 
expository  chapters  in  recent  books  (Goda,  1985,  Chapter  9;  Muga,  1984;  and  Sarpkaya  and 
Isaacson,  1981).  The  reader  is  referred  to  the  above  for  an  entry  into  the  literature.  The 
comments  here  will  be  directed  to  more  recent  results  or  to  topics  not  generally  covered  in  the 
cited  references. 

The  usual  choice  for  the  probability  law  for  individual  wave  heights  is  the  Rayleigh 
distribution,  or  some  closely  related  generalization,  such  as  the  Weibull  distribution.  If  the 
significant  wave  height  parameter,  Hs,  in  the  Rayleigh  distribution  is  determined  from  the 
traditional  4.004  times  the  area  under  the  wave  spectral  density  (Borgman,  1982a,  Equation  19), 
then  empirical  data  on  wave  heights  usually  show  a  deviation  from  the  Rayleigh  probabilities  for 
large  waves  (Forristall,  1978;  Krogstad,  1985).  An  appropriate  modification  of  the  Rayleigh 
formula  that  agrees  well  with  data  is  a  form  of  the  Weibull  distribution: 


FH(h)  =  P[Hsh]  .  l-exp[-(4h/H,)‘/P]  (175> 
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If  a  =  2  and  0  =  8,  then  the  above  formula  reduces  to  the  usual  Rayleigh  distribution  function. 
Forristall  found  that  a  =  2.126  and  0  =  8.42  were  good  choices  for  his  data.  Krogstad  reports 
values  of  a  ranging  from  2.37  to  2.50  and  values  of  0  from  12.5  to  15.6. 

An  alternative  to  the  Weibull  modification  is  suggested  by  Goda  (1978),  who  found  that 
the  significant  wave  height  is  better  estimated  by  using  3.8  times  the  area  under  the  spectral 
density  instead  of  the  narrow-band  theory  value  of  4.004.  Similarly,  Chen  (1979,  p.  19)  found 
the  bes*  multiplier  to  be  3.88  for  hurricane  Carla  and  3.82  for  hurricane  Camille  (both  Gulf  of 
Mexico  storms).  Horikawa  (1988,  p.  46)  discusses  some  of  these  results.  With  these  reduced 
values  of  Hs,  the  Rayleigh  distribution  gives  fair  results  without  modification.  Good  results  with 
the  Rayleigh  formula  arc  also  obtained  when  the  mean-square  (zero-crossing)  wave  heights  can 
be  obtained  directly  from  water  level  elevation  time  series  (Borgman,  1973b,  Figuie  7). 

The  conditional  probability  law  for  individual  wave  periods,  given  the  wave  height,  was 
shown  by  Longuet-Higgins  (1975)  to  be  normally  distributed  for  waves  with  a  narrow-band 
spectral  density.  This  result  was  shown  by  Chen  (1979)  ind  Chen,  Borgman,  and  Yfantis  (1979) 
to  hold  approximately  for  hurricane  waves  for  wave  heights  larger  than  the  significant  wave 
height.  In  particular,  Chen  found  that  (T|H  =  h)  was  approximately  normal  with  the  mean 
equal  to  0.85Tp  and  the  standard  deviation  equal  to  0. 15TpHs/h.  Here,  Tp  is  the  period 
associated  with  the  frequency  at  the  peak  or  mode  of  the  wave  spectral  density  and  Hs  is  the 
significant  wave  height.  Chen’s  results  were  based  on  an  analysis  of  wave  records  in  hurricane 
Carla  and  Camille  in  the  Gulf  of  Mexico.  Krogstad  (1985)  found  similar  results  for  North  Sea 
waves  although  the  standard  deviations  were  somewhat  larger.  Krogstad  suggests  that  this  may 
be  due  to  the  presence  of  swell  in  his  data.  Goda  (1978)  found  quite  similar  results  in  a  suite 
of  data  that  he  had  collected. 

The  joint  probability  law  of  height  and  direction  of  travel  for  individual  waves  was 
developed  by  Borgman  (1981)  and  later  extended  to  treat  the  joint  probability  law  for  height, 
period,  and  direction  (Ogbi,  1983).  The  probability  law  for  height  and  direction  appears  to  agree 
well  in  the  data  comparisons  that  have  been  made,  while  the  height-period-direction  trivariate  law 
is  much  more  speculative  and  needs  to  be  checked  more  thoroughly  against  data. 

Probability  laws  for  wave  properties  over  a  period  of  time,  such  as  the  maximum  wave 
height,  have  received  substantial  attention  in  the  literature.  In  the  multiyear  or  decade  frame  of 
reference,  probabilities  for  maximum  wave  height  become  intimately  interrelated  to  the  storm 
climatology  of  the  region.  A  very  detailed  study  of  such  topics  for  the  Gulf  of  Mexico  was 
sponsored  by  a  "consortium"  of  oil  companies  (Ward,  Borgman,  and  Cardone,  1979).  Other 
discussions  of  the  probability  law  for  long-term  maximum  wave  height  are  given  by  Ochi  (1982, 
p.  283ff.)  and  Naess  (1984).  The  probabilities  for  maximum  wave  height  over  a  time  interval 
within  which  the  random  seas  are  statistically  stationary  may  be  treated  with  extensions  from  the 
Rayleigh  distribution  (Borgman,  1973b;  Ochi,  1982,  p.  287ff.;  Goda,  1985,  p.  227)  or  from 
stochastic  process  theory  (Boccotti,  1985;  Naess,  1985). 

The  study  of  the  long-term  probabilities  for  spectral  shape  parameters  is  very  data 
intensive  and  specific  for  each  location.  Usually,  this  kind  of  study  involves  a  cooperative  effort 
between  meteorologists  for  the  wave  hindcasts,  engineers  for  relevant  structural  concerns,  and 
statisticians  for  the  probabilistic  data  analysis  (Ward,  Borgman,  and  Cardone,  1979;  Borgman 
and  Resio,  1982), 
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Wave  Group  Simulation 

The  statistics  of  wave  groups  is  really  a  topic  in  its  own  right,  with  rather  different 
mathematical  techniques  and  concepts  (Goda,  1985,  p.  230ff.),  and  will  not  be  treated  here. 
Instead,  several  comments  will  be  introduced  regarding  ways  that  wave  groups  can  be  embedded 
into  a  conditional  simulation.  If  a  short  time  history  of  wave  group  water  level  elevations  is 
available  from  measurements,  the  rest  of  the  record  can  be  conditionally  simulated  on 
reproducing  that  interval  of  data  within  the  simulation.  Thus,  one  can  produce  a  number  of 
simulations,  each  containing  exactly  the  specified  group. 

Another  way  to  force  the  occurrence  of  a  wave  group,  within  exactly  specifying  the  time 
history,  is  to  introduce  a  short  sequence  of  large  water  level  elevations  that  are  separated 
approximately  by  the  period  of  the  waves  in  the  group.  That  is,  each  large  value  is  placed  in 
the  time  record  at  a  separation  of,  say,  10  seconds  from  the  previous  and  subsequent  large 
values.  The  intermediate  time  series  values  and  the  rest  of  the  wave  record  are  then 
conditionally  simulated,  given  the  assigned  magnitudes  of  the  large  values.  The  simulated  wave 
amplitudes  in  the  group  will  be  at  least  as  large  as  the  large  values,  but  some  of  the  simulated 
intermediate  values  may  rise  a  little  height.  This  forces  the  period  of  the  waves  in  the  simulated 
group  to  vary  sightly  from  the  interval  of  separation  of  the  imposed  large  values  and  allows  the 
wave  heights  in  the  group  to  vary  somewhat  from  simulation  to  simulation. 
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Integer  in  Equation  149 
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Dimensionless  x  position  in  Equation  164 
Uniform  random  number  in  Equation  83 
Real  part  of  Amj 

Dimensionless  y  position  in  Equation  164 

x  component  of  water  particle  velocity 

y  component  of  water  particle  velocity 

z  component  of  water  particle  velocity 

Conditional  random  vector  V2,  given  that  random  vector  Vj  is 
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Negative  of  imaginary  part  of  Amj 
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Equation  28 

General  frequency  sequence  in  definition  of  the  discrete  Fourier 
transform  in  Equation  27 
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Variable  in  integration  change-of-variable  in  Equations  47  and  48 
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Delta-stretched  z  value  in  Equation  74 

Stretched  z  value  in  Equation  72 
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=  1  if  the  waves  are  traveling  toward  direction  6 

=  - 1  if  waves  are  coming  from  direction  0 
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Time  sequence  in  Theorem  C 

Discrete  Fourier  transform  of  yn  in  Theorem  C 

Complete  gamma  function 

Effective  width  in  Equation  51 
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Frequency  increment  (see  Equations  31  and  32) 
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Angle  increment 

Expectations  in  Theorem  B 
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Water  density 

Wave  property  at  zl  =  0  in  Equation  71 

Derivative  of  wave  property  with  respect  to  z  at  z  =  0  in  Equation  71 
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Variance  of  sea  surface  in  Equation  46 
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=  arctan(Y/X)  in  Theorem  A 
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Appendix 


THEOREMS 


THEOREM  A 


Let  (X,  Y)  be  random  variables  with  mean  zero  and  variance  a2.  Define 


R  =  (X2+Y2)1/2 


S>  =  arctan(Y/X) 


Then  X  and  Y  are  independent  and  normally  distributed  if  and  only  if  R  and  <p  are  independent. 
R  is  Rayleigh  distributed,  and  <p  is  uniformly  distributed  on  (0,  2k). 

Note:  The  Rayleigh  distribution  here  has  probability  density 


w  = 


(r/oz),  exp(-rz/2  o  ),  rzO 

0,  for  r  <  0 


Proof.  Let  X  and  Y  be  N(0,  a2)  and  independent. 


fx.y(x»y)  =  (1/2  it  o2 )  exp[  -  (x2  +  y2 )  /2  o2  ] 


X  =  R  cos  $ 
Y  =  R  sin® 


R-(*!+y!)ln 


A-l 


♦  —  arctan(Y/X) 


Then  the  Jacobian  of  the  transformation  is 


j  _  a(*»y> 
d(r,$) 


cos  4>  -r  sin  4> 
sin4>  rcos4» 


r 


Hence, 

f*,#(r,<W  -  (l/2it  o2)exp(-r2/2o2XJ)  -  (1/2  *)(r/2  o2)exp(-r2/2  o2) 


define  for  0  <  0  <  2rr  and  r  >  0.  But  the  density  for  a  uniform  random  variable  on  (0,  2ir) 
is: 


M4» 


1/2 11 ,  for0<<}><2n 

0,  otherwise 


So 

=  f#(4>)fR(r) 


It  follows  that  R  and  *  are  independent.  *  is  uniform  on  (0,  2ir),  and  R  is  Rayleigh 
distributed. 


THEOREM  B  (Borgman  1973a) 

N-l 

Let  i  =  fA  and  bn  =  Bmexp(12  n  mn/N)  for  n  =  0,1,2,.  ..,N-1,  be  a  covariance- 

m>0 

stationary,  mean  zero,  sequence  of  real-valued  random  variables  that  are  periodic  with  period  N 
and  a  realization  of  a  time  series  sampled  with  time  increment  At.  That  is, 

E[bJ  =  0 
and 


A-2 


C»(k)  =  Efb.b^J 

does  not  depend  on  n.  Then  if 

Af  =  !/(N  At) 

N-l 

SM(m)  =  At  Y,  C^OOapt-ilxmn/N) 

k-0 

and 

B  =  U  -  iVm 

ts  ns  in 

It  follows  that: 

I-  vo  =  VN/2  =  0: 

2.  for  0  <  m  <  N/2,  Um  =  UN.m  and  Vm  =  -VN_m; 

3.  E[Um]  =  E[VJ  =  0: 

4.  Cbb(N-k)  =  Cbb(k)  and  SBB(N-m)  =  SBB(m)  are  real-valued; 

5.  {U0,  Uj,  V,,  U2,  V2,...,U(N/2H,  V(N/2).,,  UN/2}  are  uncorrelated; 

6.  Var[U0]  =  SBB(0)Af,  Var[UN/2]  =SBB(N/2)Af,  Var[Um]  =  V-[Vm]  = 
SBBAf/2,  [for  m  =  l,2,...,(N/2)-l],  and  SBB(m)Af  =  E[|Bm|2]  for  all  m. 

Proof.  By  the  discrete  Fourier  transform  inverse, 

n  - 1 

Bm  =  T7  E  bnexp(-i2nmn/N) 
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Because  bn  is  real-valued. 

t  N1 

B*  =  T7  12  bBexp(i2«mn/N)exp(-i2niiN/N) 

N  „.0 

=  £  bnexp[-i2  7i(N-m)n/N] 

N  q.o 

~  ^N-m 

This  proves  items  1  and  2  when  combined  with  the  periodicity.  Also,  from  the  FFT  inverse, 

1  N'1 

E[Bm]  =  —  52  E[bJexp(-i2wmn/N)  =  0 
N  n.o 

This  proves  item  3.  Because, 


N-k)  =  E[bnbntN.k]  =  E[babaJ 

by  periodicity 

=  E[bn_kbJ  =  E[bnbn+k] 

by  covariance  stationarity 

=  q*® 

From  the  symmetry  on  Cbb(k), 

N-l 

SBB(m)  =  At  52  Cbb(k)exp(2nkm/N) 

k-0 


which  is  real-valued.  Thus,  item  4  is  proven. 

For  0  <  m  <  N/2  and  0  <  m’2  <  N/2,  consider: 

«i  =  E[B;Ba.]  =  {E[UmUI11.]  +E[VmVm,|  +i{E[VmUm,]  ~E[UmVm,]} 
«2  =  EtB»BJ  =  |E[U!nU|n.]+E[VmVm.}+i{E[VmUB.]-E[UmVm.]} 
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These  two  expectations  can  also  be  expressed  in  terms  of  bn  as 


V 

,  N-l  N-l 

(  \ 

m-  n-  -  mn 

/ 

=  Trf  E  E  EI*>,ba.]exp 

— i  2  7c 

/N 

& 

N2  b.o  «..0 

m'n*  +  mn, 

/ 

Now  let  k  =  n  -n,  E[bnbn.]  =  C^n’-n)  -  C^k),  and  using  the  periodicity  after  setting  n’  = 
k+n. 


\ 


w 


H-l 


N-l 


— r  52  ^(k)  exp(  -i  2  n  m*  k/N)  £  «P 

N  k-0 


n -0 


-i2  rc 


I’-m  / 
n/N 
v  +  m  ' 


But 


N-l 


52  exP 

n-0 


(  \ 

nr-m  / 

-i2  n 

U/N 

m’  +  m  j  ' 

0  ,  otherwise 


With  the  stated  constraints  on  m  and  m’,  m  =  m’  is  impossible.  Hence,  introducing  the 
definition  for  Sm, 


/ 


ei' 

««, 


(  Ssato)  Af  ) 


if  m  =  m\  m  *  0,  m  *  N/2 


if  m  *  nr 


Af\ 

BB  ,  if  m  =  nr  =  0  or  m  =  m1  =  N/2 

(SBB(m)AfJ 
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If  m  =  m\  m  #  0,  and  m  *  N/2. 


E[ui]-B[vi]  -  S^dnJAf 
E[U«V.]-Ep.V.]  -  0 

E[ui]-E[vi]  -  0 

E[U»V„]*E[U.V.]  -  0 
It  follows  that,  for  0  <  m  <  N/2, 

E[ui]  =  E[vi]  =  ymJAf 12 
Bp.VJ  -  0 

For  m  -  0  or  m  =  N/2,  Vm  =  0  by  item  1.  So  e[Uq]  =  S^fO)  Af  and  e[u^]  =  S^/O)  Af . 
Ifmf  m’, 

EpMUln.]  +  E[V1BVIn.]  =  0 
E[Vm  Um>]  -  E[Um  Vm.]  =  0 
E[UmUB.]-E[VmVm.]  =  0 
E[VmUm.]  +  E[UmVm.]  =  0 
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Therefore, 


EpmUm,]  -  Biy.VJ  =  E[YBUm.]  =  E[UBVn.]  =  0 


Finally,  if  m  =  m\ 


E[|B„f]  -  B[U**Vi]  -  S„(m)Af 


This  completes  the  proof  of  items  5  and  6. 


THEOREM  C  (Borgman,  1973a;  see  Goodman,  1957  for  continuous  version). 

Let  the  two  sequences 

N-l 

bn  =  53  Bmexp(i2nmn/N) 

m-0 

Yn  =  £  rmexp(i2nmn/N) 

B*0 

for  n  =  0,1, 2,..., N-l,  be  covariance  -  and  cross-covariance  -  stationary,  have  mean  zero,  be 
real-valued,  be  periodic  with  period  N,  and  be  a  time  series  realization  sampled  at  the  time 
increment  At.  Define  Af  =  l/(NAt)  and 


C»(k)  =  E[bnbtt+k] 
Cyy  (k)  =  E[YBYn+k] 
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C^OO  -  E[b„y..k] 
Clt(k)  -  E[T>b,.kj 


SBB(m)  SBr(m) 
S,.B(m)  Srr(m) 


N-l 


A*  E 

k«0 


C^fk) 
CYb(k)  CYY(k) 


exp(-i2  x  mk/N) 


S^fm)  -  c^fm)  -  i  q^Cm) 

B»  -  u«“iv« 
r*  = 


It  follows  that: 

1.  Both  sequences  have  the  individual  properties  in  theorem  B: 

2.  CbT00  =  CTb(N-k)  are  real-valued; 

3.  SBr(m)  =  SBr(N  -m)  =  (N  - m)  are  complex-valued; 

4.  for  0  <  m  <  N/2  and  0  <  m’  <  N/2  with  m  *  m\  (Um,  Vm)  are  uncorrelated  with 

(*m»  ♦m’) 

5.  if  0  <  m  <  N/2. 


covariance  matrix  of 


Um 

^(m) 

0 

CBrfcn)  qBr(m> 

Vm 

0 

SssCm) 

-qnr<m)  °BAm) 

K 

aas 

CBr(m) 

-qBr(m> 

© 

1 

J 

.V 

W®) 

°Br 

0  S*  j 

Af 

2 
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if  m  =  0  or  m  =  N/2, 


[um1  [Wm>  ^(“>1 

covariance  matrix  of  —  Af 

*«j  [cBr^)  Srrt®) 

S„(m)4f  -  E[B;r,] 

Proof.  Item  1  is  obviously  true.  Because, 

CYb(N-k)  =  E[Ynbn+N_k]  by  periodicity 

-  ®P>«-k  Yb] 

=  E[bn  Yn.k]  by  covariance -ctationarity 

-  C„Y(k) 

This  proves  item  2  since  the  real- value  property  follows  directly  from  the  definition. 
Continuing, 

N-l 

S^r(N-m)  =  At  £  Cby(k)exp[+i2w(N -m)k/N] 

k-0 

N-l 

=  At  E  00«P(-12*mk/N) 

k-0 

=  SBr(m) 

N-l 

=  At  52  CYb(N-k)exp[+i2nm(N-k)/N] 

k-0 
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Let  j  =  N-k  and  use  the  periodicity  CYb(N).  Then, 


H-l 

S^r(N-m)  -  At  £  C  k(j) exp[i2  x  mj/N] 

k-0 

■  SpB(m) 


This  proves  item  3. 

The  proof  of  item  4,  5  and  6  is  very  parallel  to  the  corresponding  proof  in  Theorem  B 
of  item  5  and  6.  Let 


«,  -  E[B;r„.] 

-  Ep„r..] 


and  perform  the  parallel  algebra.  This  gives: 


|  Sgr(m)  Af  j 


if  m  =  mf,  m  *  0,  m  *  N/2 


if  m  *  nr 


(sBr(m)Af> 


ifm=m’=0orm=m'  =  N/2 


The  proof  of  item  4  and  5  follows  from  the  same  steps  used  in  Theorem  B. 
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THEOREM  D 


For  X  =  1,2,  define: 


N-l  i 

b™  =  E  QLkj  exP(i  2  n  m  n/N) 

■«o  j  *  1 


bS1  -  E\X/  -  •2>*t*2) 

j-i 


A*j  “  UaJ-iVmj 


off  - 


where  {(Umj,  Vmj),  m  =  0,l,2...,N/2;  j  =  1,2,...,J}  are  independent  random  variables  with 


Voj  -  V,  h  0 


A*j  "  An- 


•J 


Var(Ump  - 


Smj  Af  A0  ,  if  m  =  0  or  m  =  N/2 

SmjAfA8/2,  if  U  m  i  N/2 


Var(VaJ)  =  SmjAfA6/2 


Here,  b*X)  is  interpreted  as  being  realization  of  a  time  series  sampled  with  time  increment  At, 
and  Af  is  defined  as: 


Af  *  1/(N  At) 
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Then  the  b*X)  are  real-valued,  and  the  covariance  relations  for  b®,  b®,  Umj,  and  VMj  are  as 
follows: 

For  X  =  1,  2: 

1.  Cov[b£X), b<A>]  =  Af  A0  £  Y  smj | Off  f  exp[ -  i2  n  m(n*  -n)/N] 

1  J  m-0  j-i 

2.  Cov  [b®, b®]  =  Af  A0  52  Y  Smj Q£j  *  Q-j  e*p[ "*  2  *  m(n> -®)/N] 

L  J  n»-0  J.l 

If  0  <  m  <  N/2: 

3.  Cov^.b™]  -  SmJ  Af A0 Re[Q^ exp(i 2 it m n/N)] 

4.  CoV[Vmj,b™]  -  SmJ Af A0 Im[Q® exp (2 * m n/N)] 

If  m  =  0  or  m  =  N/2: 

5.  Cov[Um)>b™]  -  2 S_, Af  AB  RefQ™  «p(i  2  *  m  n/N)  ] 

6.  If  0  <  m  <  N/2: 

C°v[U-l’*-l  "  -Cov[V«)-’T“]  -  S~iRe[QlV] Af  A0/2 

Cov[Vml,*™]  =  Cov[Umj,T“>]  -  SmJIm[Q«,]Af49/2 

Var[<*']  =  V*[t2>]  =  £  | On,1/ 1 2  Af  A0/2 

J"t 

Var[Umj]  =  Var[Vmj]  -  SBjAfA0/2 
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Cov[^.T«]  -  0 

7.  If  m  =  0  or  m  =  N/2,  must  be  real-valued  and: 

CovK*-]  =  S.jQ.jAfAS 

Varpmj]  -  Smj  Af  A6 

Var[<X)]  =  £  sJo^/fAfAe 

1  J  j-i  L  1 

8.  If  0  <  m  <  N/2: 

-  Cov[t™,T®]  =  E  SBj Rc[q^']* Q®] Af A0/2 

j-1 

Cov^®.?®]  -  -Cov[*®,*®]  =  E  S.jImfQS'-QSjAf  A9/2 

9.  If  m  =  0  or  m  =  N/2,  Q^}  is  real-valued  and 

Cov[®2\*®]  =  £  S.jQ^Q^AfAe 
j*» 

Proof.  The  proof  just  involves  a  careful  application  of  the  previous  theorems. 
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DISTRIBUTION  QUESTIONNAIRE 

The  Naval  Civil  Enginaering  Laboratory  is  revising  its  primary  distribution  lists. 


SUBJECT  CATEGORIES 

1  SHORE  FACILITIES 

1 A  Construction  methods  and  materials  (including  corrosion 
control,  coatings) 

1 B  Waterfront  structures  (maintenance/deterioration  control) 

1C  Utilities  (including  power  conditioning) 

1 0  Explosives  safety 
IE  Aviation  Engineering  Test  Facilities 
1 F  Fire  prevention  and  control 
1G  Antenna  technology 

1 H  Structural  analysis  and  design  (including  numerical  and 
computer  techniques) 

1 J  Protective  construction  (including  hardened  shelters,  shock 
and  vibration  studies) 

IK  Soil/rock  mechanics 
1 L  Airfields  and  pavements 
1 M  Physical  security 

2  ADVANCED  BASE  AND  AMPHIBIOUS  FACILITIES 

2A  Base  facilities  (including  shelters,  power  generation,  water 
supplies) 

2B  Expedient  roads/airfields/bridges 
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